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1.0 INTRODUCTION 


Finite-difference (FD) and finite-volume (FV) methods are very powerful 
techniques for obtaining solutions to partial differential equations that govern fluid flow 
problems. However, in order to use these methods, it is necessary to replace the spatial 
domain of the problem being studied by a finite number of discrete points known as 
grid points. The process of replacing a spatial domain by a system of grid points is 
referred to as grid generation. Grid generation is a very important part of FD and FV 
methods because the system of grid points used strongly affects the accuracy, 
efficiency, and ease with which these methods generate solutions. In some instances, 
the ability or inability to generate an “acceptable” grid system determines whether FD 


or FV methods can or cannot be used. 



Even though tremendous advances have been made in grid generation 
techniques during the past fifteen years (refs. 1 to 10), the generation of acceptable 
grid systems for geometrically complex three-dimensional spatial domains remains a 
difficult problem. Recently, a very efficient and versatile computer program, called 
GRID2D/3D, has been developed which can generate grid systems inside complex- 
shaped two- and three-dimensional (2- and 3-D) spatial domains. GRID2D/3D is so 
efficient that it is configured to run on PCs or PC compatible computers, though it can 
also be used on workstations and mainframes. The high efficiency of GRID2D/3D 
makes it especially useful for spatial domains that deform with time. This is because 
for such spatial domains, a different grid system must be generated at each time level, 
and the number of time levels can be thousands or more. This technical memorandum 
describes the theory and method behind GRID2D/3D and the types of grid systems 
that it can generate. Part 2 of this technical memorandum (to be published under a 
separate cover) will contain the program, GRID2D/3D, and a user’s manual. 

1.1 Types of Grid Systems That Can Be Generated by GRID2D/3D 

Eiseman and Erlebacher (ref. 9) classified all possible grid systems that can be 
used by FD and FV methods as follows. At the broadest level, a grid system can be 
classified as structured, unstructured, or mixed depending upon how the grid points 
are connected to each other (fig. 1-1). A structured grid system, in turn, can be 
classified as a single grid or a composite grid. A single grid is one that is based on a 
single boundary-fitted coordinate system, whereas a composite grid is made up of two 
or more single grids patched together with each single grid having a different 
boundary-fitted coordinate system. Depending upon how the different single grids are 
patched together, a composite grid can further be classified as completely 
discontinuous, partially discontinuous, partially continuous, or completely continuous 
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(fig. 1-2). The continuity or discontinuity referred to here is concerned with that of 
the different boundary-fitted coordinate systems at locations where they are patched 
together in a composite grid. 

Of the grid systems mentioned above, the unstructured grid system is the most 
versatile and the easiest to generate, especially for complicated-shaped spatial domains. 
But, the use of unstructured grid systems with FD and FV methods is still at a state 
of development (refs. 11 to 16). Presently, FD and FV methods almost exclusively use 
structured grid systems, and that is the type of grid system GRID2D/3D generates. 

When a structured grid system is used with a FD or a FV method to obtain 
solutions to fluid flow problems, the structured grid system generated by GRID2D/3D 
or any other computer program should satisfy a number of conditions and they are 

1. The total number of grid points in the grid system should be kept to the 
minimum needed for the FD or FV method to yield solutions of the desired 
accuracy. This condition is important for computational efficiency and can 
be achieved by clustering grid points in regions where they are needed (e.g., 
regions where gradients of the flow are large) and scattering them elsewhere. 

2. One set of grid lines (coordinate lines of the boundary-fitted coordinate 
system) always should coincide with the boundary of the spatial domain 
regardless of the geometric complexity or motion of that boundary (i.e., the 
grid system should be boundary conforming). This condition is important 
because it enables FD and FV methods to implement boundary conditions 
easily and accurately for geometrically complex and/or deforming spatial 
domains. 

3. Grid lines that intersect a boundary should intersect that boundary 
perpendicularly so that derivative boundary conditions can be implemented 
more easily and accurately. At the interior of the spatial domain, the angle 
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of intersection between grid lines only needs to be nearly orthogonal (i.e., 
between 45 and 135 degrees). 

4. The spacings between grid points should change slowly from a region where 
grid points are concentrated to a region where grid points are sparsely 
distributed, especially in regions where gradients of the flow are large. This 
condition is important because Fourier components which make up the 
solution reflect and refract at interfaces where grid spacings change. 

5. One set of grid lines should align with the flow direction. This condition is 
important for convection dominated flows when the aspect ratio of the 
control volume about each grid point is very high and/or when the thin layer 
Navier-Stokes equations are used to study such flows. 

For complicated 2- and 3-D flows within geometrically complex spatial domains, 
the flow field varies considerably from one region to another. Thus, it is usually not 
possible to generate a single grid that would satisfy all of the above conditions at every 
part of the spatial domain. For such fluid flow problems, it is often necessary to 
generate a number of different single grids, each of which satisfies the above five 
conditions at a different part of the spatial domain. These single grids are then 
patched together to form a composite grid. 

As noted earlier, depending upon how the different single grids are patched 
together, a composite grid can be completely discontinuous, partially discontinuous, 
partially continuous or completely continuous. In general, the more discontinuous a 
composite grid is, the easier it is to generate that grid and the harder it is to use that 
grid to obtain solutions. 

Since the computer program, GRID2D/3D, is capable of generating single grids 
as well as the different types of composite grids, we briefly discuss below the 
advantages and disadvantages of the various types of composite grids in order to know 
when a specific type should be used. 
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1.1.1 Completely Discontinuous Composite Grids. - The major advantage of 
completely discontinuous composite grids, such as the chimera grid (refs. 17 to 21 and 
fig. l-2(a)), is that they are the easiest to generate. To illustrate how chimera grids 
are generated, consider the spatial domain for the flow past an entire aircraft. To 
generate a chimera grid for such a spatial domain, all one has to do is generate a series 
of single grids, one about each component of the aircraft; for example one about the 
fuselage, another about the wing, still another about the nacelle, and so on. The 
patching process simply involves laying each single grid over the appropriate 
component of the aircraft, deciding the amount of overlap of different single grids, and 
ensuring that the entire spatial domain is filled with grid points. Since the geometry 
for each single grid can be made relatively simple and patching is trivial, the grid 
generation process is straightforward. 

Another important advantage of this type of grid is that the structure of each 
single grid can be different from each other; for example, one single grid may have a C- 
C structure, while another may have an 0-0 or an O-H structure. Thus, it is possible 
to optimize each single grid for a different part of the spatial domain. 

Still another important advantage of this type of grid is that it is the easiest to 
do local grid refinement. For a chimera grid, one can refine the grid at any location by 
simply generating a very fine single grid and then overlaying it wherever desired. 

Also, this type of grid can easily be applied to problems in which one or more 
objects are moving relative to another object, such as the launching of missiles from an 
aircraft (ref. 20). For such problems, the completely discontinuous composite grid may 
be the best type of grid to use. 

Finally, since composite grids are composed of a series of single grids, it is 
possible to do computations on one single grid at a time. This will reduce computer 
memory requirements considerably since only information on one single grid needs to 
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reside in the computer at any one time. This advantage is shared by all composite 
grids, continuous or discontinuous. 

The major disadvantage of completely discontinuous composite grids is that it 
is more difficult to obtain solutions on such grids when using FD and FV methods than 
the other types of composite grids. This is because interpolation and averaging 
schemes are needed to transfer information from one single grid to another when and 
wherever two or more single grids overlap (refs. 18 to 21). Also, the schemes must be 
developed to ensure that boundary conditions are implemented correctly for the entire 
problem and that certain properties, such as the conservative and the transportive 
properties, are maintained in regions where two or more single grids overlap. 

1.1.2 Partially Discontinuous Composite Grids. — Partially discontinuous 
composite grids (fig. l-2(b)) are generated in the following manner. First, the spatial 
domain of the problem being studied is partitioned into a number of nonoverlapping, 
contiguous zones or blocks. Next, a single grid is generated within each zone. Finally, 
patching of the single grids simply involves putting each of the single grids into its 
respective zone. 

Thus, partially discontinuous composite (PDC) grids differ from completely 
discontinuous composite (CDC) grids in that the single grids of PDC grids do not 
overlap each other. However, PDC and CDC grids have two important similarities. 
First, each single grid in both cases can have a structure that is different from each 
other. Second, the number of grid points in each single grid can be different from 
each other. Because of these two important similarities, the major advantages of PDC 
grids are very similar to those of the CDC grids. However, since single grids in a PDC 
grid do not overlap each other, PDC grids are somewhat more difficult to generate but 
are easier to use than CDC grids. 
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Examples of PDC grids include the “zonal” or “patched” grids (refs. 22 to 
25). For such grids, the main difficulty is the implementation of boundary conditions 
at the interfaces where the different single grids of the PDC grid meet. References 22 
to 25 describe a procedure for implementing such boundary conditions in a way that 
would ensure maintenance of the conservative property. 

1.1.3 Partially and Completely Continuous Composite Grids. — Completely 
continuous composite (CCC) grids are grid systems in which all grid lines (i.e., 
coordinate lines of the boundary-fitted coordinate system) and all of their derivatives of 
every order are continuous at all interfaces where different single grids meet. In 
general, it is not necessary to construct CCC grids. Typically, FD and FV methods 
only require continuity of the grid lines and their first and, occasionally, second-order 
derivatives at the interfaces where different single grids meet. Composite grid systems 
with this limited degree of continuity are referred to as partially continuous composite 
(PCC) grids. 

Figure l-2(c) shows a PCC grid in which grid lines are all continuous, but first- 
order derivatives of the grid lines have discontinuities. It can readily be seen in that 
figure that the slope of the grid lines and the spacing between the grid lines change 
suddenly at the interface where the two single grids meet. Figure l-2(d) shows a PCC 
grid in which the grid lines and their first-order derivatives are continuous everywhere 
including the interface where the two single grids meet. Such PCC grids have the same 
appearance as CCC grids. 

The major advantage of PCC grids of the type shown in figure l-2(d) is that 
this is the easiest grid system for FD and FV methods to use. This is because 
boundary conditions can be implemented easily at interfaces where different single grids 
meet. In fact, it is not even necessary to treat the interfaces where different single 
grids meet as boundaries since computations can be carried across them. For PCC 
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grids, the complete spatial domain of the problem can be mapped onto a single 
transformed domain, even though different boundary-fitted coordinate systems have 
been used in different parts of the spatial domain (fig. 1-3). 

Here, it is important to note that not all grid systems which appear to be 
continuous are continuous. Figure 1-4 shows a composite grid that appears to be 
continuous but belongs to the PDC grids because it is impossible to map the entire 
spatial domain onto one transformed domain. 

The major disadvantage of PCC grids is that they are the most difficult 
to generate when compared to CDC and PDC grids. Another disadvantage of PCC 
grids is that the structure and number of grid points in each single grid must satisfy 
certain compatibility conditions to ensure continuity. These compatibility conditions 
make it more difficult to optimize each single grid for a specific area. It also makes it 
more difficult to do local grid refinement. 

Thus, there are many spatial domains for which it is extremely difficult, if not 
impossible, to generate a PCC grid that is acceptable. However, when it is possible, 
then the generation of such grids is worthwhile because of the ease with which they 
can be used. References 26-31 show a number of examples of how to construct PCC 
grids for complex-shaped 2- and 3-D spatial domains. 

1.2 Grid Generation Methods Used in GRID2D/3D and Why 

In the previous section, we discussed the various types of grid systems that can 
be generated by GRID2D/3D and when a specific type should be used. In this 
section, we briefly outline the different types of advanced grid generation techniques 
and give reasons why GRID2D/3D is based on one class of methods. 

All grid generation techniques can be divided into two major classes — 
differential equation methods and algebraic methods. Differential equation methods 
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generate grid systems by solving a system of partial differential equations (PDEs) 
which describes how grid points are to be distributed within the spatial domain. 
Examples of differential equation methods include methods based on elliptic PDEs 
(e.g., harmonic mapping, 2- and 3-D quasilinear systems), methods based on hyperbolic 
PDEs (e.g., orthogonal trajectories and field methods), and methods based on parabolic 
PDEs (refs. 2, 3, and 8). Most of these methods require a significant amount of 
computational effort since the systems of PDEs that must be solved are quasilinear and 
often as complicated as the PDEs that govern the fluid flow problem. This is 
especially true when using these methods to generate grid systems in 3-D spatial 
domains and in spatial domains that deform with time. 

Algebraic methods generate grid systems by interpolating between boundaries 
of the spatial domain. Since no PDE needs to be solved in the grid generation process, 
algebraic grid generation methods are computationally much more efficient than 
differential equation methods. Examples of algebraic grid generation methods include 
shearing transformations (ref. 32) and transfinite interpolation methods (refs. 33 to 
35). Transfinite interpolation methods include the Two-Boundary Method (refs. 36 
and 37), Four-Boundary Method (refs. 38 and 39), Six-Boundary Method (refs. 39 and 
40), and multisurface methods (refs. 41 to 43). 

Whether one uses a differential equation method or an algebraic method, the 
grid generation process is always iterative. This is because the “acceptable” grid 
system is arrived at via trial and error after generating a series of unsatisfactory grids. 
The effort of the iterative process is, of course, compounded many times for spatial 
domains which can deform with time since, for such domains, a different grid system is 
needed for each time level and the number of time levels can be thousands or more. 
Hence, the efficiency of the grid generation process is extremely important for problems 
with 3-D spatial domains and for problems in which the spatial domain can deform. 
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Since GRID2D/3D is intended for complex-shaped 2- and 3-D spatial domains 
and for spatial domains that can deform with time, GRID2D/3D generates grid 
systems by using algebraic grid generation methods. Depending upon the complexity 
and dimensionality of the spatial domain, GRID2D/3D uses one of the following three 
methods, all of which are very similar and are based on transfinite interpolation: the 
Two-Boundary Method, the Four-Boundary Method, and the Six-Boundary Method. 
These methods were chosen because of their high efficiency and their ability to provide 
very precise controls over the distribution of grid points in the spatial domain when 
used in conjunction with stretching functions (refs. 44 and 45). Also, these methods 
can generate grid lines that intersect boundaries orthogonally. 

By using these algebraic grid generation methods, GRID2D/3D generates single 
grids with grid lines that are continuous and differentiable everywhere up to the 
second-order. GRID2D/3D generates composite grids by patching together two or 
more single grids. The patching can be discontinuous or continuous. For continuous 
composite grids, the grid lines are continuous and differentiable everywhere up to the 
second-order except at interfaces where different single grids meet. At interfaces where 
different single grids meet, the grid lines are only differentiable up to the first-order. 

In order to use the Two-, Four-, and Six-Boundary Methods to generate grid 
systems, the boundaries of the spatial domains must be represented mathematically in 
parametric form. This is a difficult problem for complicated-shaped spatial domains 
because the boundaries of such domains are complicated as well. In GRID2D/3D, 
parametric equations for boundary curves of 2-D spatial domains are generated by 
either spline interpolation (ref. 46) or tension-spline interpolation (ref. 47). Parametric 
equations for boundary surfaces of 3-D spatial domains can be generated by a number 
of techniques including linear Coons’ interpolation (ref. 33), bidirectional spline 
interpolation (refs. 48 and 49), and bihyperbolic spline interpolation (ref. 50). In 
GRID2D/3D, parametric equations for these 3-D surfaces are generated by either 
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linear Coon’s interpolation, bihyperbolic spline interpolation, or a new technique 
referred to as 3-D bidirectional Hermite interpolation. 

The details of the algebraic grid generation methods used in GRID2D/3D are 
given in Sections 2.0 and 3.0. Examples of single and composite grids generated by 
GRID2D/3D are given in Section 4.0. A listing of the computer program, 

GRID2D/3D, and a user’s manual will be found in Part 2 of this technical 
memorandum which will be published under a separate cover. 
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2.0 THE TWO-. FOUR-, AND SIX-BOUNDARY METHODS OF ALGEBRAIC GRID GENERATION 


As noted in Section 1.0, GRID2D/3D generates grid systems by using one of the 
following three algebraic grid generation methods: the Two-Boundary Method, the Four- 
Boundary Method, and the Six-Boundary Method. All three of these methods are based on a 
general technique known as transfinite interpolation (refs. 33 and 35). The Two- Boundary 
Method, described by Smith (ref. 36) and Yang and Shih (ref. 37), is intended for problems in 
which it is only necessary to map correctly two arbitrary-shaped boundaries of the spatial 
domain. For problems in which it is necessary to map correctly all of the boundaries of the 
spatial domain or at least more than two of them, then the Four-Boundary Method or the Six- 
Boundary Method, described by Vinokur and Lombard (ref. 38), Rizzi and Eriksson (ref. 39), 
and Eriksson (ref. 40), can be used. All three of these methods can generate grid systems in 3- 
D spatial domains; the Two- and Four- Boundary Methods can also generate grid systems in 2- 
D spatial domains. 

In this section, the details of these methods are described by using each of them to 
generate either a 2-D or a 3-D grid system. Our step-by-step descriptions of the methods 
follow closely those of Yang and Shih (ref. 37) and Shih (ref. 51). It is our intention to present 
the methods in a clear manner so that the reader might easily implement any one of the three 
methods. 

2.1 The Two- Boundary Method 

As mentioned previously, the Two-Boundary Method is intended for problems in which 
it is only necessary to map correctly two arbitrary-shaped boundaries of the spatial domain. 


12 



There are a number of problems for which this holds (Section 6-4-3 of ref. 51); a few example 
problems with such spatial domains are shown in figure 2-1. 

Here, we note that the Two-Boundary Method can correctly map all of the boundaries 
of a spatial domain if the remaining boundaries are straight lines in the 2-D case or flat 
surfaces in the 3-D case. Figures 2- 1(a) and 2-l(d) illustrate these cases. 

The Two- Boundary Method involves the following eight major steps (refs. 37 and 51): 

1. Define the nature of the coordinate transformation. 

2. Select a time-stretching function. 

3. Select the two boundaries of the spatial domain that must be mapped 
correctly. These two boundaries cannot touch each other at any point. 

4. Describe the two boundaries selected in Step 3 in parametric form. 

5. Define curves that connect the two boundaries using transfinite 
interpolation. 

6. Discretize the domain (i.e., replace the continuous domain of the problem by 
time levels and grid points). 

7. Control the distribution of the grid points with stretching functions. 

8. Calculate the metric coefficients needed by the FD or the FV method to 
obtain solutions. 

The details of these eight steps of the Two-Boundary Method are described below by 
generating a grid system in a 3-D, deforming spatial domain shown in figure 2-2(a) for the 
problem of compressible flow through a converging-diverging channel. The spatial domain of 
interest is the region bounded by surfaces 1 through 6 in figure 2-2(a). The spatial domain 
deforms because surfaces 1 and 2 deform with time. 

Step 1 - Define the Coordinate Transformation. The first step of the Two- 
Boundary Method is to define the coordinate transformation between the coordinate system of 
the spatial domain and the boundary-fitted coordinate system of the transformed domain. For 
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3-D spatial domains in which grid points are allowed to move, grid generation involves the 
determination of the following coordinate transformation: 

(x,y,z,t) *-> (£,y,C,r) (2.1a) 

or, more specifically, 

(2.1b) 
(2.1c) 
(2- Id) 
(2-le) 

where x, y, z , and t represent the coordinate system of the spatial domain and rj y £, and r 
represent the boundary-fitted coordinate system of some transformed domain (fig. 2-2). 

Step 2 - Select a Time-Stretching Function. The next step is to define a 
relationship between t and r. For our example, we set t equal to r; that is, 

t = r (2.2) 

Thus, no time-stretching function is used. Time stretching may be useful when variable time- 
step sizes are used with FD or FV schemes that involve information at more than two time 
levels. 

Step 3 - Select Two Boundaries of the Spatial Domain. The third step is to 
select the two boundaries of the spatial domain that are to be mapped correctly. These two 
boundaries must not intersect each other at any point. For the spatial domain of figure 2- 
2(a), we select boundary surfaces 1 and 2. 


t ~ t{r) 
x = x(t,r}X,T) 

y = 
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Since £, 77 , £ and r represent a boundary-fitted coordinate system, boundary surfaces of 
the spatial domain in the x-y-z-t coordinate system must correspond to coordinate planes in 
the £-r)-(,-r coordinate system. We choose surfaces 1 and 2 to correspond to coordinate planes 
77 = 0 and 77 = 1 , respectively (fig. 2 - 2 ); that is, 


Xj = 4Z, *7=0, C, T) = X 2 (£,C,r) 

(2.3a) 

Y 1 = y(£, rj=0, C, r) = Y 2 (£,<,r) 

(2.3b) 

Z 2 = *?=0, C, r) = Z 2 U,C,r) 

(2.3c) 

x 2 = *7 = 1, C» T ) = X 2 (f,C,r) 

(2.4a) 

Y 2 = y(£, y=l, C, r ) = Y 2 ({,(,r) 

(2.4b) 

z 2 = z(Z, 77 = 1 , C, r) = Z 2 (££,t) 

(2.4c) 


Here, X 2 , Y 2 , and Z 2 are the r-, y- , and ^-coordinates of surface 1, and X 2 , Y 2 , and Z 2 
are the x-, y-, and z-coordinates of surface 2. The remaining four boundaries — surfaces 3, 4, 5, 
and 6 -- are mapped to coordinate planes £ = 0 , £ = 1 , £ = 0 , and £ = 1 , respectively (fig. 
2-2). 

Step 4 - Describe the Two Boundaries Selected in Parametric Form. Once the two 
boundaries have been selected, the next step is to represent these two boundaries in parametric 
form as suggested by the form of equations (2.3) and (2.4). Equations (2.3) and (2.4) also tell 
us that the three parameters which must be used to describe surfaces 1 and 2 are £, £, and r. 
Keep in mind that had we chosen different coordinate planes in the transformed domain for 
our boundaries, these parameters could have been different. 

In this example, we assume that the parametric equations describing surfaces 1 and 2 
are given and they are 

x, = ( / L x 

Yj — A sin(wr) [1 — cos(2 7 r £)] [1 — cos(2 7 r £)] 

Z i = C / L y 
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(2.5a) 

(2.5b) 

(2.5c) 



and 

x 2 = ( / L x 

Y 2 = Ly — A sin(wr) [1 — cos(2 tt £)] [1 — cos(2 tt C)j 

Z 2 = C / Ly 


(2.6a) 

(2.6b) 

(2.6c) 


where A, w, L x , Ly ? and L z are given constants. 

For most problems, parametric equations describing the boundaries of spatial domains 
are not given. What is usually given are sets of coordinates which describe the positions of a 
finite number of discrete points located on the boundary, and numerical methods must be used 
to generate the required parametric equations. The numerical methods used in GRID2D/3D 
for this purpose are described in Section 3.0. 

Step 5 - Define Curves Connecting Boundaries Using Transfinite Interpolation. A 
number of different transfinite interpolation techniques can be used to derive curves which 
connect the two boundaries chosen in Step 3. Here, we consider two such methods: transfinite 
interpolation based on Lagrange interpolation and transfinite interpolation based on Hermite 
interpolation. 

Transfinite interpolation based on Lagrange interpolation is also known as linearly 
blended transfinite interpolation. When this technique is used to generate connecting curves 
between surfaces 1 and 2, the resulting curves have the following functional form: 

= X 2 (£,C,r) l^rj) + X 2 (£,C,t") l 2 (*}) (2.7a) 

y{Z,V,(,r) = Y J (^,C,r) l^rj) + Y 5 (£,C,r) l 2 (rf) (2.7b) 

= Z j(U,t) h(v) + Z 2 (U,r) l 2 (r}) (2.7c) 

where X^, Yj, and Z 2 given by equation (2.5) describe surface 1, and X^, Y^, and Z 2 given by 
equation (2.6) describe surface 2. The functions / 2 and l 2 are blending functions, and, for the 
linear case, they have the functional forjn of first degree Lagrange interpolating polynomials. 
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As such, they connect two points — one on each surface having the same £, C, and r values — 
and are constrained by the following expressions: 

hi 7 ) = o) = i hi 7 ) = 1) = o 

I 2 (n = 0) = 0 i 2 (v = 1) = 1 

With these constraints, / 1 and l 2 become 

hi 7 )) = 1 ~ V 
hi 7 )) = V 

Substitution of equation (2.8) into equation (2.7) yields the desired linear connecting curves 

= XjK.C.t) (1 — i?) + X 2 («,r) n (2.9a) 

!«>•»,<, r) = Y 1 ((,(,r) (1-r,) + Y 2 «,<,r) 9 (2-9b) 

= Zjtf.C.r) (1-r;) + Z 2 ((,(,r) , (2.9c) 

It is often desirable for connecting curves to intersect boundaries orthogonally so that 
derivative boundary conditions can be implemented accurately. Since the curves described by 
equation (2.9) are straight lines, they will not, in general, intersect boundaries orthogonally. 
One way to remedy this is to use transfinite interpolation based on Hermite interpolation to 
form the connecting curves. Transfinite interpolation based on Hermite interpolation allows 
specification of the derivatives at the end points of the curves, enabling one to force 
orthogonality at the boundary. When this method is used to generate connecting curves 
between surfaces 1 and 2, the resulting cubic curves have the following functional form: 


(2.8a) 

(2.8b) 
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= X^.C.r) h^rj) + X 2 (£,C,r) h 2 (tj) 


dx(£, 77 = 0 , C,r) 5x(£,i?=l,C,r) , ( , 

+ ^ W) + ^ h M 

y(t>*hC,T) = Y faCr) h^rj) + Y 2 (S£,t) h 2 (r}) 

, 0y(£,*?=O,C,r) k / V <9y(£,7?=l,<,r) 

+ 5^ h sW + 5^ M>» 

z{Z,yX,r) = ZjiZX ,r) h^rj) + Z p (£,< ,t) h 2 (rj ) 

, dz(^r]=OX,r) u dz(£,r)= l,C,r) ( , 

+ 3T M?) + Jfn h M 


(2.10a) 


(2.10b) 


(2.10c) 


where X ; , Y ; , and Z i given by equation (2.5) describe surface 1, and X 2 , Y 2 , and Z 2 given by 
equation (2.6) describe surface 2. The functions hj, h 2 , h 3 , and h 4 are blending functions 
which connect two points — one on each surface having the same £, and r values. They are 
constrained by the following expressions: 


hj{v = 0) 

= 1 

h^T) = 1) 

- 0 

dh^j] = 0) 
dr) 

- 0 

dhj(r] = 1) 
dr) 

= 0 

h 2 (v = 0) 

= 0 

h 2 (v = 1) 

= 1 

dh 2 (r} — 0) 

= 0 

dh 2 (r] = 1) 

= 0 

dr) 

drj 

M 7 ? = °) 

= 0 

h 3 (v = !) 

= 0 

dh 3 {t] = 0) 

_ 1 

dh 3 (t] = 1) 

— 0 

dr) 

— X 

dr) 


h 4 (n = o) 

= 0 

h 4 (r) = 1) 

= 0 

dh 4 (, = 0) 

= 0 

dh 4 (rj = 1) 

— 1 

dr) 

dr) 



With these constraints, h p h 2 , /i 3 , and h 4 become (ref. 36) 


h^Tj) — 2rj 3 — 3i) 2 + 1 
h 2 (r)) = ~2 V 3 + Sr) 2 
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(2.11a) 

(2.11b) 



(2.11c) 
(2. lid) 


h 3 (v) = T } 3 - 2r] 2 + T] 
M*?) = V 3 ~ V 2 


We choose the values for 


dx(Z,Ti=0£,T) dx(Z,r)=l,C,T) dy(t>V = 0X,r) 


and 


dr) ’ dr) ’ dr) 

()y(£ i f T ) 

- so that the connecting curves given by equation (2.10) will intersect surfaces 1 

and 2 orthogonally. For surface 1, this will occur when the cross product of n (a vector 
normal to surface 1) and e n (the vector tangent to the connecting curve) is zero. This will be 
the case when 


dx(C,r)~ 0,C,r) _ 


dY jdZ 2 

d ZflY^ 

dr) 

■ dC dC 

dC dt) 

dy(C,y=0X,r) _ 

-/<-,({, (,r)( 

dX 1 dZ 1 


dr) 

dC dt 

dC dt ) 

II 

■O 

o 

II ^ 

^(ec,r)( 

dX 1 dY 1 
dC dt 

dY 2 dX 2 \ 
dC dt) 


(2.12a) 


(2.12b) 


(2.12c) 


Similarly, the connecting curves will intersect surface 2 orthogonally when 

(2.13a) 
(2.13b) 
(2.13c) 

K 1 (tX> T ) and ^(£>C r ) m equations (2.12) and (2.13) are known as the “K factors” 
and are chosen by trial and error so that no overlapping of the connecting curves takes place 
in the interior of the spatial domain. For our problem, the “K factors” are constants given by 


dx(£,r) = l,C,r) 
dr) 


I< 


>(£>C, r )( 


dY 2 dZ 2 dZ 2 dY 2 


dt d C dt d( 


) 


dy(^,r)-lX,r) _ /dX 2 dZ 2 dZ 2 dX 2 \ 

drj ~ SC d£ dC ) 


dz(C,r) = l,C,r) _ 
dr) 


I< 




dX 2 dY 2 dY 2 dX 2 
dt dC dt d( 


) 
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= K 2 (U,t) = 0.2 


n-\L , ,dXi(U,r) dY 2 (U,T) dX 2 (U,r) A dY 2 {Us) . . /n 

The values of L -^ , L -^ , , and in equations (2.12) 

and (2.13) can easily be found by analytically differentiating equations (2.5) and (2.6) or by 

using finite-difference formulas. Thus, substitution of equations (2.11), (2.12), and (2.13) into 

equation (2.10) yields the desired cubic connecting curves based on Hermite interpolation. 


Step 6 - Discretize the Domain. Steps 1 through 5 above describe how we map the 
x-y-z-t coordinate system onto the £-t)-(-t coordinate system. Having done this, we now need 
to discretize the domain in the £-r)-(-r coordinate system; that is, we need to replace the 
continuous domain by time levels and grid points. 

For our problem, we replace the time domain by equally incremented time levels; that 
is, 


r n = n A r , n = 0, 1, 2, ... (2.14) 

where r n denotes the time at time level n, and A r denotes the constant time-step size. 

We also replace the spatial domain in the coordinate system by IL X JL x KL 

equally spaced grid points (fig. 2-3(b)). The locations of these grid points are given by the 
ordered triples where 


*< = 

= (*— 1) Af , 

i — 1> 2, 

..., IL 

(2.15a) 


: (j-1) At) , 

3 = 1, 2, 

..., JL 

(2.15b) 

£ k = 

= (*-1) AC , 

k — 1, 2, 

..., I<L 

(2.15c) 


_ 1 
“ IL- 1 

^ = jrh 

1 “ KL- 1 

(2.15d) 


If we substitute equation (2.15) into equation (2.10), we can obtain the locations of the 
grid points in the x-y-z-t coordinate system. Figure 2-3(a) shows the system of grid points for 
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our problem in the x-y-z-t coordinate system obtained by using equations (2.10) and (2.15) at 
one value of r. 

Step 7 - Control the Distribution of Grid Points. At this point, we need to 
examine the grid system shown in figure 2-3(a) and ask, “Is the distribution of grid points 
satisfactory?” In order to answer this question, we need to consider the physics of the problem 
for which the grid is generated. For accurate solutions, grid points should be clustered in 
regions of the spatial domain where sharp gradients in the dependent variables exist. Such 
clustering can be achieved by the use of stretching functions (refs. 44 and 45). 

For compressible flow within the 3-D spatial domain shown in figure 2-2(a), we expect 
steep gradients near the walls (surfaces 1, 2, 3, and 4). Grid points can be clustered near 
surfaces 1 and 2 by replacing r) in equation (2.10) by the following stretching function 
expression: 


(Pv + 1, 

!((/?„+l)/(/?„-l)) ( 

2, " 1) — /S.i + 1 

2 { 

! + (/? >? + !)/ (Pii — 1) 


\ 


(2.16) 


where /3 n is a constant greater than unity. More clustering takes place in the rj direction near 
y — O and T] = 1 as j3 v approaches unity. 

In a similar manner, grid points can be clustered near surfaces 3 and 4 by replacing £ 
in equation (2.10) by 


(^+i)((^+i)/(^-i)) (2e ~ 1) -/? f +i 

2{i+(/? e +i)/(/? { -i)) ( ^' 1) } 


(2.17) 


where (3 ^ is a constant greater than unity that acts in the £ direction as does in the rj 
direction. 
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The new distribution of grid points in the x-y-z-t coordinate system after stretching is 


shown in figure 2-4. 

Step 8 - Calculate Metric Coefficients. Once we obtain a satisfactory 

distribution of grid points in the x-y-z-t coordinate system, we are ready to calculate the metric 
coefficients which are needed to obtain finite-difference solutions to the partial differential 
equations governing the problem. The metric coefficients appear in the governing equations 
when they are transformed from the coordinate system of the spatial domain ( x-y-z-t in our 
example) to the boundary-fitted coordinate system of the transformed domain (£- 77 - (-t). With 
the coordinate transformation described by equation ( 2 . 1 ), the metric coefficients which appear 
are r t , £ t , r) t , ( t , £ x , Vx, Cx, £y, V y, Cy> £*, V z, and <*. These metric coefficients can be 
calculated using the following expressions (refs. 2, 3, and 8 ): 


r t = 1 (2.18a) 

tx = (y v z c -y c z n )/ J (2.18b) 

£ y = -{xyz^-x^zy)/ J (2.18c) 

iz = (x^y^-x^yr,)/ J (2.18d) 

it = -[ x T (Vn*<-ycZt,) ~ y^xyz^-x^) + z T (x„t/ ( -x ( y„) )/ J (2.18e) 

r\x = -(y^ c -y c z f )/7 (2.18f) 

Vy = (x^-x^)/ J (2.18g) 

r]z = -(x f y c -a; c y e )/ J (2.18h) 

V t = [ * T (ytZ ( -y c z € ) - y r (x ( z c -x <Z ^ + z T (x^y ( -x c y f ) ]/ J (2.18i) 

Cx - {y^-yyZ^)/ J ( 2 . 1 8 j ) 

Cy = -(x^z v -x n z ( )/ J (2.18k) 

Cz = (x ( yr,-x n y^/ J (2.181) 

Ct = - [ x r(y^ z n~yn z 0 - y T ( x t z rt- x n z t) + (x^Vn-x^y^) ]/ j (2.18m) 


where J is the Jacobian given by 
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J = x^y^-y^zr,) - + x ( (y f *„-y„^) 


(2.18n) 


The derivative terms x v , x^ y y,,, y^, and in equation (2.18) can be evaluated 

either analytically by differentiating equation (2.10) or numerically by using finite-difference 
formulas. The correct way to evaluate the metric coefficients depends on how the governing 
equations written in the boundary-fitted coordinate system are cast. If the governing 
equations are cast in strong conservation-law form, then the metric coefficient must be 
evaluated numerically and the finite-difference formulas used to evaluate them must be the 
same as those used in the finite-difference method of solution. If the governing equations are 
cast in weak conservation-law form, then no conditions are imposed on the method for 
evaluating the metric coefficients themselves but there is a condition on how the derivatives of 
metric coefficients can be evaluated. Finally, if the chain-rule conservation-law form is used, 
then no conditions are imposed on how metric coefficients and their derivatives are evaluated. 
This important topic is addressed in references 51 to 54. 

2.2 The Four-Boundary Method 


The Four-Boundary Method for generating grid points is intended for situations where 
four boundaries of a spatial domain need to be mapped correctly from the spatial domain to 
the transformed domain. The Four-Boundary Method (refs. 38, 39, 40, and 51) is an 

extension of the Two-Boundary Method described in the previous section, and, as such, 
consists of the same eight major steps with minor variations, namely: 

1. Define the nature of the coordinate transformation. 

2. Select a time-stretching function. 

3. Select the four boundaries of the spatial domain that must be mapped 
correctly. 
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4. Describe the four boundaries selected in Step 3 in parametric form. 

5. Map the spatial domain to a transformed domain. 

6 . Discretize the domain (i.e., replace the continuous domain of the problem 
with time levels and grid points). 

7. Control the distribution of the grid points with stretching functions. 

8 . Calculate the metric coefficients needed by the FD or the FV method to 
obtain solutions. 

We will describe each of the eight steps of the Four-Boundary Method by generating a 
grid system inside the spatial domain shown in figure 2-5(a). This spatial domain is the cross- 
section of a steel bar whose temperature distribution we wish to determine. The bar has a 
uniform cross-section along its length, we shall consider only the 2-D region bounded by curves 
1 through 4 (fig. 2-5(a)). 

Step 1 - Define the Coordinate Transformation. For the 2-D, non-deforming 

spatial domain shown in figure 2 - 5 (a), we seek a coordinate transformation of the form 

(x,y,t) <-► (£,T 7 ,r) (2.19a) 

or, more specifically, 

t = t(r) (2.19b) 

x — ^£, 17 ) (2.19c) 

y = y(Z,v) (2.i9d) 

Here, x, y, and t comprise the coordinate system of the spatial domain, and £, rj, and r 
comprise the boundary-fitted coordinate system of the transformed domain (figs. 2-5(a) and 
2-5(b)). 
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Step 2 - Select a Time-Stretching Function. A time-stretching function is not used 
and t is set equal to r as shown by equation (2.2). 

Step 3 - Select Four Boundaries of the Spatial Domain. Since the spatial 
domain of figure 2-5(a) has only four boundaries, all four boundaries are selected. We choose 
curves 1, 2, 3, and 4 to correspond to coordinate lines 77 — 0, 17 = 1, £ = 0, and £ = 1, respectively 
(fig. 2-5); that is, 


x ; = *<<;, ’?=«) = Xj(«) ( 2 - 20a ) 

Yj = y(«, y=0) = Y,({) (2.20b) 

X 2 = x(£, y = l) = X 2 (0 (2.20c) 

Y 2 = y({, 0 = 1) = Y 2 (4) (2.20d) 

X 5 = z«=0, 0) = X 2 (o) (2.20e) 

Y 3 = y(£=0, 0 ) = Y 2 (o) (2.20f) 

X, = *(« = 1, 0 ) = X*(ij) (2.20g) 

Y 4 = y(£ = l> 0) = Y 4 (o) (2.20h) 


Here, X ; and are the x- and ^-coordinates of curve 1; X 2 and Y 2 are the x- and y- 
coordinates of curve 2; X 3 and Y 3 are the x- and y-coordinates of curve 3; and X 4 and Y 4 
are the x- and y-coordinates of curve 4. 

Step 4 - Describe the Four Boundaries Selected in Parametric Form. Having selected 
four boundaries, we now need to represent these boundaries in parametric form as suggested 
by equation (2.20). For this problem, information about the four boundary curves is given in 
the form of a set of discrete points which lie along the curves. Thus, interpolation techniques 
must be used to generate the parametric equations of the boundary curves. In GRID2D/3D, 
either spline or tension spline interpolation can be used. Here, tension spline described in 
Section 3.0 is used to obtain parametric equations for the four curves in terms of the 
parameter £ for curves 1 and 2 and in terms of t) for curves 3 and 4. Figure 2-6 shows the 
curve approximations thus generated. 
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Step 5 - Map the Spatial Domain. The Four-Boundary Method maps the spatial 
domain to the transformed domain in two steps. The first step is essentially the same as Step 
5 of the Two-Boundary Method in which we select two of the four boundaries which do not 
touch each other and which have been described parametrically using the same parameter 
(either f or ij). As in Step 5 of the Two-Boundary Method, curves that connect these two 
boundaries are specified by using Hermite transfinite interpolation. When this step is 
completed, the two boundaries that were selected will be mapped correctly, but the other two 
boundaries will in general be mapped incorrectly. To remedy this, a second step is performed 
where the mapping constructed during the first step is modified so that the other two 
boundaries will also be mapped correctly. 

In our example, we will first define curves that connect curves 1 and 2 such that only 
curves 1 and 2 will be mapped correctly. Afterwards, we will modify the connecting curves so 
that curves 3 and 4 will also be mapped correctly. 

Curves which connect curves 1 and 2 are described by the following Hermite 
interpolation expressions: 


At*) = X 2 (0 hfo) + X 2 (0 h 2 (rj) 


+ 


ds(£,q=0) 

dr) 


h 3 (r)) + 


dx(£,7?=l) 

dr) 


M 7 ?) 


(2.21a) 


= Y 2 (0 h t ( V ) Y,(0 h 2 ( V ) 


+ 


dy{j,r)~ 0 ) 

dr) 


h 3 {r}) + 


dy{£,,r)~ 1 ) 

drj 


h 4 to 


(2.21b) 


Here, h v h 2 , h 3 , and h 4 are given by equation (2.11). The values of the partial derivatives in 
equation (2.21) will be given shortly. 
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By using equation (2.21), we map curves 1 and 2 in the x-y-t coordinate system to 

coordinate lines r)—0 and 77 = 1 in the £-r)-r coordinate system (fig. 2-7). Note, however, that 

curves 3 and 4 have not been mapped to coordinate lines f = 0 and £ = 1. Instead, connecting 
curves 3' and 4' have been generated in the spatial domain, and it is these curves which have 
been mapped to £ = 0 and £ = 1 in the transformed domain. In order to map curves 3 and 4 in 
the x-y-t coordinate system to coordinate lines f=0 and £ = 1 in the £-r)-T coordinate system, 
we must adjust curves 3' and 4' so that they coincide with curves 3 and 4. With this in mind, 
we define two quantities — Ax and Ay — such that 

x(t,T}) = x'((,Ti) + A x(Z,T]) (2.22a) 

y(Z,v) — y'(Z,v) + Ay(£,y) ( 2 . 22 b) 

will map curves 3 and 4 to £=0 and £=1, respectively. Here, x> and y' are given by equation 

( 2 . 21 ). 

Ax and Ay are given by the following Hermite interpolation expressions: 


A x(£,y) = [x(£ = 0,y) - x , (^ = 0,y)] h 5 (£) 

+ [x(£ — 1,77) - x , (^ = l,y)] /i 6 (£) 

. r dz(£ = 0 , 77 ) ^^(^=0,77) 

+ l ^ J 

, r dx{^~l ,77) < 9 x'(£ = 1,t/) 1 (e , 

+ [ — ^ J h 8tt) 

Ay(C,v) = [y(^— 0,77) - y'(f= 0 , 77 )] h 5 (£) 

+ [y(£ = Uv) - y'(f = M)] h 6 (Z) 

rdy(^-0,T)) dy'(Z=0,r}) 

+ l ^ J 

dy(S = l,rj) dy'(Z = l,T}) 

+ [ — dl dl — 1 h *W 


(2.23a) 


(2.23b) 
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where 


<M£= 0,r/) L . , dz((=0,ri=0) , , . fl*(£=0,ij=l) 

MO + MO 


a 2 ^=o,^=o) , a 2 x((= o,„=i) 

+ h M> atari + h4{r,) Sidr, 

(2.23c) 

dx'(i = 1,0 , / \ dx(i = l,Tj=Q) , , , dx(i 1 , 1 , 1) 

d( - MO at + 2W di 


(n) d 2 4i= 1 . 1 = 0 ) + h (n) d 2 4i =M=i) 

+ WJ d { dr) + n 4W) didf) 

(2.23d) 

dy'(Z=0,y) L / \ dy(£ = 0,r?=0) , , , dy(t = 0,T) = l) 

at - mO d( + m ,} di 


a 2 y((=o,n=0) d 2 y(i=o,n=i) 

+ h M dtdrj + h M> dtdrj 

(2.23e) 

dy'(Z= M) , f x 0y(£ = l,»/=O) , x = M=l) 


1 h rn^ 2y (^ = 1,77 = °) 1 A M ^y(£ = M = l) 

(2.23f) 

h 5 (0 = 2^ - 3i 2 + 1 

(2.24a) 

MO = -2f 3 + 3£ 2 

(2.24b) 

MO = i 3 -2i 2 + ( 

(2.24c) 

MO = t 3 -t 2 

(2.24d) 


If we substitute equations (2.23) and (2.24) into equation (2.22), then we obtain the 
desired expressions for x(£,t/) and y(£ ? f/) which describe the mapping between the spatial 
domain and the transformed domain. 

We still need to specify the derivative terms in equations (2.23) and (2.24). Similar to 
the Two- Boundary Method, the first-order derivative terms are chosen so that the connecting 
curves will intersect the boundaries orthogonally. This time, however, the spatial domain is 
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two dimensional so that we must use the dot product instead of the cross product to specify 
orthogonality at a boundary. Connecting curves will intersect curve 1 orthogonally when the 
dot product of e ( (a vector tangent to curve 1) and (the vector tangent to the connecting 
curve) is zero. It can be shown that this will be the case when 


dxjZ,T]=0) _ 
drj 


-*i(0 


0 Yj(O 

dZ 


and 


dy{Z,r}~ 0) 
dr} 


*i(0 


flXj(Q 

dZ 


(2.25) 


Following this line of reasoning, the expressions for the first-order derivative terms in 
equations (2.23) and (2.24) are given below: 


dx(£,7? = 0) _ 
dr} 


dy{Z,v= 0) 

’ dr} 


(2.26a) 

dx{Z,r} = l) _ 
dr\ 

-K»r a f 

dy(Z,r} = l) _ 
' dr} 

KM “f 

(2.26b) 

dx{Z~0,r)) 

dZ 


dy(Z-0,y) 

-W”£" 

(2.26c) 

dx(Z = l,r}) 
dZ 


dy(Z = l,r}) 

' dZ 


(2.26d) 


For our example , Ki(Z) an d were chosen to be equal to 0.3, while K 3 (r}) and K 4 (r}) were 

chosen to be equal to 0.1. 

Methods for determining the second-order derivative terms present in equation (2.24) 
are given in reference 40. In our example, these terms are all set equal to zero. 

Step 6 - Discretize the Domain. We discretize the domain in the coordinate 

system by replacing the temporal domain with equally incremented time levels and by 
replacing the spatial domain with IL x JL equally spaced grid points (fig. 2-8(b)). 
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The time levels are described by equation (2.14). The grid points are located at (£,-,»?•) 


where 


€<= (*-i) A* , 

i= 1, 2, ..., IL 

(2.27a) 

Vj= 0-1) At? , j= 1, 2, ..., JL 

(2.27b) 

= IL—l 

^ = JL- 1 

(2.27c) 


If we substitute equations (2.23), (2.24), (2.26), and (2.27) into equation (2.22), then 
we can obtain the locations of the grid points in the x-y-t coordinate system. Figure 2-8(a) 
shows the system of grid points for our problem in the x-y-t coordinate system obtained by 
using equation (2.22). 

Step 7 - Control the Distribution of Grid Points. In this example, we choose not 
to use stretching functions to redistribute the grid points within the spatial domain. If 
redistribution is desired, then the procedure described in the previous example can be followed. 

Step 8 - Calculate Metric Coefficients. For our 2-D, non-deforming spatial domain, 
the metric coefficients which need to be evaluated are r t , £*:, r } x , and r]y. These metric 

coefficients can be evaluated by using the following equations: 

T t = 1 
£x = Vri/ j 

Vx = -y^/J 
€y = —x v / J 

T)y - Xf/J 

where J is again the Jacobian but this time is given by 

J = x^y v - x n y ( (2.28f) 


(2.28a) 

(2.28b) 

(2.28c) 

(2.28d) 

(2.28e) 
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As before, the partial derivative terms in equation (2.28) can be evaluated either 
analytically or numerically by using finite-difference formulas depending upon how the 
governing equations written in the boundary-fitted coordinate system are cast. 

2.3 The Six-Boundary Method 

The Six-Boundary Method for generating grid points is intended for 3-D spatial 
domains in which six boundaries of the spatial domain need to be mapped correctly from the 
spatial domain to the transformed domain. The Six-Boundary Method (refs. 39, 40, and 51) is 
an extension of the Two-Boundary Method and the Four-Boundary Method described in the 
previous two sections, and, as such, consists of the same eight major steps with minor 
variations, namely: 

1. Define the nature of the coordinate transformation. 

2. Select a time-stretching function. 

3. Select the six boundaries of the spatial domain that must be mapped 
correctly. 

4. Describe the six boundaries selected in Step 3 in parametric form. 

5. Map the spatial domain to a transformed domain. 

6. Discretize the domain (i.e., replace the continuous domain of the problem 
with time levels and grid points). 

7. Control the distribution of the grid points with stretching functions. 

8. Calculate the metric coefficients needed by the FD or the FV method to 
obtain solutions. 
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We shall describe each of the eight steps of the Six-Boundary Method by generating a 
grid system within the deforming, spatial domain shown in figure 2-9(a). 

Step 1 - Define the Coordinate Transformation. For the 3-D, deforming spatial 

domain shown in figure 2-9(a), we seek a coordinate transformation of the form 


(x,y,z,f) «-► (£,r?,C,r) 


(2-29a) 


or, more specifically, 


t = i(r) 

(2-29b) 

x - x(£, 77 ,C,r) 

(2.29c) 

y = y(Z,r),(,T) 

(2.29d) 

z = >r) 

(2.29e) 


Here, x, y, z, and t comprise the coordinate system of the spatial domain, and £, 77 , £, and t 
make up the boundary-fitted coordinate system of the transformed domain (figs. 2 - 9 (a) and 2 - 
9(b)). 

Step 2 - Select a Time-Stretching Function. A time-stretching function is not used 
and t is set equal to r as shown by equation ( 2 . 2 ). 

Step 3 - Select Six Boundaries of the Spatial Domain. Since the spatial domain 
of figure 2-9(a) has only six boundaries, all six boundaries are selected. We choose surfaces 1, 
2, 3, 4, 5, and 6 to correspond to coordinate planes 77 =0, 77 = 1, £=0, f = l, £= 0, and ( = 1, 
respectively (fig. 2-9); that is, 


X, = n =0£, t) = Xjff.C.r) (2.30a) 

Y j = y((, > 7 = 0 , C,r) = Yjtf.O) (2.30b) 

X 2 = x(t, n=U(,r) = X 2 («,r) (2.30c) 

Y 2 = y(f, >7 = 1, c,r) = Y 2 (£,C,r) (2.30d) 
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X 3 = 1(4=0, T),C,r) = X 3 (i),(,T ) (2.30e) 

Y, = y(£ — 0, r = Y 3 (t ),<, r) (2.30f) 

X„ = *({=1, r,X,r) = X 4 (^,(,r) (2.30g) 

Y< = *(4 = 1, nXs) = Y 4 (i),(,r) (2.30h) 

X 5 = a«, i),< = 0,r) = X s ((,t),T) (2.30i) 

Y s = y(A, rjX = 0,r) = Y 5 «,ij,t) (2.30j) 

X 6 = *(4, V ,C = 1,1-) = X 6 ((, V ,t) (2.30k) 

Y e = y(i, nX = l,r) = y 6 ((,V,r) (2.301) 


Here, Xj and Yj are the x- and y - coordinates of surface 1; X 2 and Y 2 are the x- and y- 
coordinates of surface 2; X 3 and are the x- and y-coordinates of surface 3; X 4 and Y 4 are 
the x- and y-coordinates of surface 4; X 5 and Y 5 are the x- and y-coordinates of surface 5; X 6 
and Y^ are the x- and y-coordinates of surface 6. 

Step 4 - Describe the Four Boundaries Selected in Parametric Form. Having selected 
six boundaries, we now represent these boundaries in parametric form. For this problem, 
information about the six boundary surfaces is given in the form of a set of discrete points 
which lie along the surfaces. Thus, interpolation techniques must be used to generate the 
parametric equations of the boundary surfaces. In GRID2D/3D, a new technique referred to 
as 3-D bidirectional Hermite Interpolation can be used (Section 3.0). 

Step 5 - Map the Spatial Domain. We map the spatial domain to a transformed 
domain in three steps. In the first step, correctly map two of the six boundaries by using the 
Two-Boundary Method. In the second step, correct the mapping completed in the first step 
by ensuring two more boundaries are mapped correctly (same as Step 5 of the Four-Boundary 
Method). Finally, in the third step, correct the mapping completed in the second step by 
ensuring the remaining two boundaries are mapped correctly. 

Surfaces 1 and 2 will be mapped correctly by the following Hermite interpolation 
expressions: 
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= X^^C.r) hjir]) + X 2 (^,C,r) h 2 {rj) 


dx(Z,r)=0£,T) u dx(t,T)=l£ ,r) 

^ MO + ^ MO 


y'(£,*?,<,r) = Yj(^,C,r) ^( 77 ) -f Y 2 U,<,t) MO 

, <9y(^,»? : =o,C,r) L , v . 0yU,*7=i><»O . / n 

+ ^ MO + ^ MO 


Here, h 2 , h 2 > h 3 , and h 4 are given by equation (2.11). 

Surfaces 1, 2, 3, and 4 will be mapped correctly by the following equations: 

x"(f,77,C,r) = 

- y'(^,r/,C,r) + Ay'(£,7?,C,0 
where x 7 and y' are given by equation (2-31) and 

As'foi^O = [X 5 0?*CO - ^U=0,t?,C,r)] MO 
+ [X*(i/,C,r) - * , (€ = l,i?,C,r)l MO 
, r 5x(^ = 0,7?,C,r) a^=0,i7,C ,0i t 

+ 5? 5? 1 MO 

, r ax(^ = l,7/,C,r) 0x'(£ = l,i 7 ,<,r) 1 t 

+ [ ^ J MO 

Ay'(£,t 7 ,C,r) = [Y 5 (rj,C,r) “ y'(£=0,7?,C,O) MO 
+ [Y 4 (vX,r) - y'(£ = M,C,r)] MO 

. r M£=0»»?>C»O ay'(^=0,»/,C,r) 1 

+ l 3£ 5? j MO 

. f dy(t=ltf,C,r) dy , (£=l i T)£,r). . 

+ [ 5| 3? J 

In the above equations, h 5 , h 6 , h 7 , and are given by equation (2.24). 


(2.31a) 


(2.31b) 


(2.32a) 

(2.32b) 


(2.33a) 


(2.33b) 
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All six surfaces of the spatial domain shown in figure 2-9(a) — namely, surfaces 1, 2, 3, 
4, 5, and 6 — will be mapped correctly by the following equations: 

00 = *"(0»7>C»0 + Ax"(£,rj£,r) (2.34a) 

y(OM» r ) = y"(0M>0 + Ay"(0M>0 (2.34b) 

where x ” and y" are given by equation (2-32) and 


Ax"(0*7 ,00 - [X 5 (Z,T},r) - z"(OM=0,O] MO 

+ PM0*7,0 - £"(0*?,C=1,0] h 10 (£) 

, r MO»7,<=0,O ^^"(^,»7,C=0,7-) 1 L 

+ l 5? j ^O 

, r M0M = l>0 ^a/ , (^,i7,C=l,T') 1 L 

+ l i M(0 


(2.35a) 


Ay'(0*?,00 = [Y 5 (y,C,r) - y'(£=0,7j,C,O] h 5 (t) 

+ [Y 4 (i 7 ,C.f) - y'(£=l,»?,00] MO 


, r <9y(£ = 0y,00 dy'(£=0,»?,<,r)i L „ x 

+ 1 M 5? ] 7 ^ } 

, r dy(f = l,77,C,r) <V(£ = M,00i L ,, x 

+ 1 dl dl J hs ^ } 

MO - 2<0 - 3C 2 + i 
M(0 = -2C 5 + 3<* 

M(0 = C 3 ~ 20 + C 
M(0 = C 5 - C 2 


(2.35b) 

(2.36a) 

(2.36b) 

(2.36c) 

(2.36d) 


If we substitute equations (2.31) to (2.33) into equation (2.34), then we obtain the 
desired expressions for £(0*7,00 and y(0*?,00 which describe the mapping between the 
spatial domain and the transformed domain. 
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We still need to specify the derivative terms in equations (2.31), (2.33), and (2.35). 
Similar to the Two- and Four-Boundary Methods, the first-order derivative terms are chosen 
so that the connecting curves will intersect the boundaries orthogonally. Methods for 
determining the second-order derivative terms are given in reference 40. In our example, these 
terms are all set equal to zero. 

Step 6 - Discretize the Domain. We discretize the domain in the ^-rj-C-r coordinate 
system by replacing the temporal domain with equally incremented time levels and by 
replacing the spatial domain with IL x JL x KL equally spaced grid points. 

The time levels are described by equation (2.14). The grid points are located at 
(tirtj’Ck) where V are S iven h y equation (2.15). 

If we substitute equation (2.14) and (2.15) into equation (2.34), we can obtain the 
locations of the grid points in the x-y-t coordinate system. 

Step 7 - Control the Distribution of Grid Points. In this example, we choose not 
to use stretching functions to redistribute the grid points within the spatial domain. If 
redistribution is desired, then the procedure described in the example in Section 2.1 can be 
followed. 

Step 8 - Calculate Metric Coefficients. For our 3-D, deforming spatial domain, the 
metric coefficients which need to be evaluated are given by equation (2.18). 

2.4 Additional Remarks 


We conclude this section by mentioning three problems which must be dealt with when 
using the Two-, Four, and Six-Boundary Methods. 

First, slope discontinuities present in the boundaries of spatial domains will propagate 
into the interior of the grid systems generated by using these methods. These discontinuities 
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in the slopes of the grid lines are undesirable since they can lead to errors in the solution. To 
correct for this, some technique should be used to smooth the grid. One way to smooth a grid 
system with slope discontinuity is to apply a Laplacian operator to the region near the 
discontinuity. This method is described in Section 4.0. 

Second, care must be taken when choosing the “K factors” and the stretching functions 
for the Two-, Four-, and Six- Boundary Methods. Large “K factors” tend to produce grid lines 
with more curvature. Such grid lines often overlap in the spatial domain or, at least, can form 
a grid system which is very skewed. In general, numerical values in the “K factors” and the 
stretching functions are arrived at in an iterative manner. First, a grid is generated by using 
one set of inputs for the “K factors” and the stretching functions. Next, that grid is plotted 
(GRID2D/3D contains a graphics program to plot grid systems that it generates) and 
inspected visually. Based on that inspection, the inputs are modified accordingly. This process 
repeats until a satisfactory grid has been obtained. Since GRID2D/3D is highly efficient, an 
acceptable grid system can be generated within a short time. 

Last, when connecting curves are discretized to form grid points, the orthogonality 
which was forced at the boundaries may be lost. Figure 2-10 illustrates this point. Between 
grid points a and b, curve c is approximated by line segment a-b. The original 90° angle, a, 
between boundary d and curve c has been replaced by angle (3 between boundary d and line 
segment a-b. In order that /I more nearly approximate a, two things can be done. First, 
stretching functions can be used to move point a closer to point 6. Second, larger “K factors” 
can be utilized to force the effect of orthogonality further into the domain along curve c. In 
this way, orthogonality between the grid lines and the boundary curves can be maintained 
after the discretization of the spatial domain. 
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3.0 METHODS FOR GENERATING PARAMETRIC 


REPRESENTATION OF BOUNDARIES 


In Section 2.0, it was shown that in order to use the Two-, Four-, and Six-Boundary 
Methods, it is necessary to represent the boundaries of the spatial domain in parametric form. 
These boundaries are curves for two-dimensional spatial domains and surfaces for three- 
dimensional ones. In this section, methods for representing curves and surfaces in parametric 
form are described. 


3.1 Parametric Representation of Curves and Surfaces 


Curves and surfaces can be described mathematically in several different ways. For 
example, a curve in the xy-plane can be represented by 

y ~ A x ) 

or 

A x ,y) - o 

Alternatively, we can describe the same curve by 

x = g(s) y = f [y(s)] - h(s ) (3-3) 

Equation (3.3) is a parametric representation of the curve in terms of a parameter, s. The 
choice of s is rather arbitrary, the only restriction being that 5 must increase monotonically 
along the curve. Here, we note that all curves and surfaces can be represented by parametric 
equations and that there is no one unique way of representing a curve or surface in parametric 
form. 


(3.1) 

(3.2) 
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In grid generation, information about a curve or surface is given either by an analytical 
expression such as equation (3.1) or by a set of coordinates which describe the locations of a 
finite number of discrete points on the curve or surface. When information about a curve or 
surface is provided in the form of an analytical expression (such as eq. (3.1)), it is a 
straightforward matter to generate a set of parametric equations (such as eq. (3.3)) for the 
curve or surface. When information about a curve or surface is given by a finite number of 
discrete points, then some type of interpolation schemes must first be used to approximate the 
curve or surface by an analytical expression before it can be represented in parametric form. 
The following subsections present methods to obtain parametric equations for curves and 
surfaces given as sets of discrete points. 


3.2 Approximation of Curves in Two and Three Dimensions 


For 2-D spatial domains, Lagrange interpolation is usually unsatisfactory for 
representing curves because it produces curves which may oscillate wildly when the number of 
discrete points (henceforth referred to as nodal points) is large (e.g., ref. 55). Hermite 
interpolation is usually impractical as well, since it requires information about derivative values 
at the nodal points which is seldom available (e.g., ref. 56). Least-squares regression is not 
computationally efficient for large sets of nodal points and yields curves which do not, in 
general, pass through each nodal point (e.g., ref. 57). Spline interpolation, on the other hand, 
yields curves that do pass through each nodal point and are well behaved (i.e., they do not 
oscillate wildly). Spline interpolation uses a different function (typically a low-degree 
polynomial) between successive nodal points (fig. 3-1). The resulting piecewise curve is made 
smooth by enforcing continuity of as many derivatives as possible at the nodal points. This 
technique can easily be used to approximate curves which are functions of two or three spatial 
coordinates. Because of the aforementioned attractive features of spline interpolation, this 


39 



method will be used here to approximate curves in both two and three dimensions. We now 
consider two types of splines — cubic splines and tension splines. 


3.2.1 Cubic Spline Interpolation. — The objective in cubic spline interpolation is to 
derive a different third-degree polynomial for each interval between successive nodal points, 
subject to the condition that the polynomials be connected piecewise to form a smooth curve 
that is continuous up to the second-order derivative. Since we require the curve to be written 
in parametric form, the equation describing this curve in an arbitrary interval, «, can be 
written as (fig. 3-1) 

X.(s) = A iS 3 + B,s 2 + C t s + D { (3.4) 

Here, X represents any spatial coordinate (e.g., £, y, or z in the Cartesian system), and s is a 
parameter which was mentioned previously and which will be discussed in more depth later in 
this subsection. A { , B { , <7-, and D { are coefficients which vary from interval to interval. For 
n+1 nodal points, there are n intervals and 4 n such coefficients; thus, 4n conditions are 
required in order to evaluate these coefficients. These conditions are summarized below: 

1. The spline curve must pass through the nodal points. This gives 2n 

conditions. 

2. The first- and second-order derivatives of the spline curve must be 

continuous at all interior nodal points (i.e., nodal points 1 through n— 1; fig. 3-1). 
This gives 2n — 2 conditions. 

3. Two conditions control the behavior of the spline curve at the two end 

nodal points. These conditions will be discussed later in this subsection. 

We could use these 4n conditions together with equation (3.4) to obtain 4n 

simultaneous equations in 4n unknowns. However, a shortcut is possible which requires the 
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solution of only n—1 simultaneous equations (e.g., refs. 55-59). Since the shortcut method is 
much more efficient, it is used in GRID2D/3D. 

We begin the shortcut method by noting that since the curve for each interval is a 
cubic, the second-order derivative within each interval is a straight line. Thus, we can write 


X "(s\ — X"(s 1 ^ S — + X"(V)— — S -~~- 


(3.5) 


where X^ ,f (s) is the value of the second-order derivative of X at any location s within the 
interval (fig. 3-1). The parameter s takes on the values s i _ 1 and s i at the beginning and end of 
the interval, respectively. 


By integrating equation (3.5) twice with respect to s, we obtain an expression for .X^s) 
which contains two constants of integration. We solve for these constants by imposing the 
conditions that X^s) must equal X^j,) at s i t and X t (s) must equal X(s i ) at s t -. This yields 




+ 


+ 


6 J (i 

■ ' 5 r ■' • 3) 


(3.6) 


The only unknowns in the above equation are the second-order derivatives at the beginning 
and end of the interval - X"^) and X"(s,). Values of these second-order derivative terms 
can be found by invoking the condition that the first derivatives at the nodal points must be 
continuous: 


xa*i) = x M) 


(3.7) 
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Equation (3.6) can be differentiated to give an expression for X/(s). If we do this for 
both the and the I th intervals and set the two results equal to each other according to 

equation (3.7), we obtain the following result: 


( s i~ s i-l)^"( s i-i) + + ( s i+i s i)X"(Si+i) 




(3.8) 


If equation (3.8) is applied at all interior nodal points, then we obtain n— 1 
simultaneous equations in n+1 unknown second-order derivative terms. This system of 
simultaneous equations has a tridiagonal coefficient matrix and requires two more conditions 
before it can be solved. Here is where our two additional conditions regarding the end nodal 
points are necessary. There are several conditions which can be imposed at the end nodal 
points. We mention only the following two: 

1. Cyclic End Conditions. If the curve being approximated is cyclic (i.e., the 
two end nodal points represent the same point in space as in a closed curve), then it 
is appropriate to write 

X'(s 0 ) = X'{s n ) and X"{s 0 ) = X”(s n ) (3.9) 

With this assumption, the coefficient matrix of the system of equations mentioned 
previously is now cyclic tridiagonal, and the system becomes n-f 1 equations in n-f-1 
unknowns. An example of a cubic spline with cyclic end conditions is shown in 
figure 3- 2(a). 

2. Natural End Conditions. If we allow the second-order derivative of the spline 
curve to be zero at the two end nodal points, then we minimize the total curvature 
of the spline curve. This is the so-called natural spline, and the end conditions are 
given by 
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X"(s 0 ) = X"(s n ) = 0 


(3.10) 


In this case, the system of equations mentioned previously reduces to n— 1 equations 
in n— 1 unknowns. An example of a cubic spline with natural end conditions is 
shown in figure 3-2(b). 

In either case, the system of equations in A"(s-) can be solved efficiently using simple 
computational algorithms (e.g., refs. 55 to 59). Once the X"(s t -)’s are known, values for X at 
any location s along the curve can be found using equation (3.6). Recall that X in equation 
(3.6) represents any spatial coordinate — x, y or z (comments under equation (3.4)). Thus, 
the required parametric equation is given by x = X($), y = T(s) for plane curves and x = 
A(s), y = K(s), z = Z(s) for twisted curves and surfaces. 

We now return to the question of how to choose the parameter, s. One simple and 
effective way is to let s represent arc length along the curve. Of course, we do not know the 
actual arc length a priori, but an adequate estimate can be obtained by summing linear 
distances between nodal points. Thus, for a two-dimensional curve in the xy-plane (i.e., a curve 
which is a function of only two spatial coordinates — say x and y), we have s 0 — 0, Sj = s 0 
+ s 2 = s, + + and so on. 


3.2.2 Tension Spline Interpolation. — Cubic spline interpolation produces curves that 
are well behaved in most cases; however, in cases where the curve to be approximated 
possesses extreme curvature, curves generated using cubic spline interpolation often have 
unwanted wiggles. This phenomenon is illustrated in figure 3-3. One solution to this problem 
is to put the spline in tension. One might visualize this as pulling on the ends of the cubic 
spline to “straighten out” any unwanted wiggles. We add tension to the cubic spline by 
replacing equation (3.5) by 
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(3.11) 


+ [*'^)-<^(s,)]g^ 


where u is the tension parameter, the range of which is given by 0 < <x < oo. The spline curve 
created by adding tension to the cubic spline is known as a tension spline (ref. 56). The 
tension spline tends toward a linear spline as a is increased and approaches a cubic spline for 
values of a near zero. 


As with the cubic spline, we integrate equation (3.11) twice with respect to s and 
require that the curve pass through the proper nodal points to give 


XM = 


X"i s t -i) sinh[a(s i — s)] 
a 2 «tn*[<r(5 1 . — 5 t ._ 2 )] 








cr 2 sin%(»r s i.i)l ' a 2 


(3.12) 


Upon demanding continuity of the first derivative at the interior nodal points, we have 




( S »“ 5 t-j) S * n %( S i“ S i-i)lj <? 2 
r o-cosh[<r(s--s-, J )] l 


+ 


+ 


+ 


\ -Sj- 1 )] ( S i~ S i-l) 

acoshltriSj+j-Sj)] _ l 

smh[<r(s.^ i -5-)] { s i+i~ s i)j <r 2 

f 1 a 

\( s i+i~ s i) sinh [<r( s i+i~ s i)}f <r 2 


' Xis^-XjSj) _ X(s i )-X{8 i , 1 ) 

( s i+l~ s i) ( s i~ s i-l) J 


(3.13) 


Applying equation (3.13) at all interior nodal points again results in a system of n— 1 
equations in n + 1 unknowns with a tridiagonal coefficient matrix. Two conditions are required 
at the end nodal points to make solution of the system possible. The cyclic and natural end 
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conditions presented for the cubic spline can be used for the tension spline as well with X"(s 
replaced by — o ,2 X(s i ). In either case, the resulting system of simultaneous linear 

equations can be solved by using simple numerical algorithms. 

Since tension splines involve hyperbolic functions, they are computationally less 
efficient than cubic splines; however, excellent results can be obtained for cases where cubic 
splines prove unsatisfactory. Figure 3-4 shows that a tension spline with cr — 10 exhibits none 
of the wiggles visually present in the cubic spline of figure 3-3. The same set of nodal points 
was used to generate the spline curves shown in figures 3-3 and 3-4. 


3.3 Approximation of Surfaces in Three Dimensions 


For three-dimensional spatial domains, the boundaries of the domain are surfaces. As 
one might expect, methods for describing surfaces based on interpolating between specified 
nodal points which lie on the surface are more complex than their two-dimensional 
counterparts. Part of the additional difficulty lies in the fact that there are several ways in 
which information about discrete points on a surface may be given. With curves, it is 
sufficient to start at one end of the curve and specify nodal points at random intervals as one 
moves along the length of the curve. This is not the case for surfaces. One cannot just 
randomly pick points on a surface and expect to efficiently construct an approximation of the 
surface using an interpolation method unless the points are chosen in an organized manner. 
Here, we consider the following two structured methods for specifying points on a surface: 

Case 1. Only points which lie along the edges of the surface are given. These edges 
exist as twisted (3-D) or plane (2-D) curves. This case is applicable to surfaces whose 
interior regions are similar to their edge contours. Examples of such surfaces are 
shown in figure 3-5. 
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Case 2. Points on the surface are given on a grid as shown in figure 3-6. In figure 3-6, 


the z-coordinates of surface H are given on a rectangular grid in the Cartesian system. 
This case is applicable to any surface but is especially well-suited to surfaces whose 
interior regions differ greatly from their edge contours. Examples of such surfaces are 
shown in figure 3-7. 

In the next three subsections, we present three methods for approximating boundary 
surfaces by interpolating between points which lie on the surfaces. Section 3.3.1 describes a 
method, referred to as transfinite interpolation with bilinear blending (also known as linear 
Coon’s interpolation), which can be used when the nodal points are given in the format of 
Case 1 above. Section 3.3.2 presents a new method, referred to as three-dimensional 
bidirectional Hermite interpolation, which can also be used when the nodal points are given in 
the format of Case 1. Section 3.3.3 presents a method, referred to as parametric bihyperbolic 
spline interpolation, which can be used when the nodal points are given in the format of Case 2 
above. 


3.3.1 Transfinite Interpolation With Bilinear Blending - To illustrate how transfinite 
interpolation with bilinear blending approximates surfaces, consider the boundary surface 
shown in figure 3-8. That boundary surface is bounded by four twisted curves. Information 
about the four twisted curves may be given in analytical form, or they may be put into 
analytical form from coordinates of the nodal points located on the curves by using either 
cubic spline or tension spline interpolation described in Section 3.2. In either case, we end up 
with parametric equations describing each of the curves in the following form: 


A 2 — X^Sj) 
Yi=Yi( S i) 
Z j — % i( s l) 


X2 — X 2 {S2) 

Y 2 =Y 2 (s 2 ) 

Y>2 — ^2^ s 2) 


X 3 =X 3 (s 3 ) 
Y 3 =Y 3 (s 3 ) 
Z3 — Z 3 (s 3 ) 
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X 4 = X 4 (s 4 ) 
Y 4 =Y 4 {s 4 ) 
Z 4 — Z 4 (s 4 ) 



Here, X- y Y^, and Z i describe curve i where i — 1, 2, 3, 4, and s- represents the approximate 
arc length along curve i as discussed in Section 3.3.1. 

Recall that for the Two-, Four- and Six-Boundary Methods, boundary surfaces in the 
spatial domain are mapped onto coordinate planes in the transformed domain. Here, we map 
the boundary surface shown in figure 3-8 onto the coordinate plane shown in figure 3-9 located 
at £ = 0 in the transformed domain. Accordingly, we relate s^, s 2 , s 3 , and s 4 to £ and r\ as 
follows: 


= 5 j (^,? 7=0,C=0) = s ; (£) (3.14a) 

= s*(£,iy=l,C=0) = s 2 (£) (3.14b) 

s 3 = Sj(£=0,i7, C=0) = s 3 (rj) (3.14c) 

s 4 = s 4 (£=1,v,C=0) - s 4 (r) i) (3.14d) 


By using equation (3.14) above, we linearly interpolate between curves 1 and 2 to 
obtain the following equations: 


x 12 ((,tiX=0) = (1 — 17 ) *,(0 + n X 2 (() (3.15a) 

Vi2((,nX=0) = (1-1?) y,(0 + v Y 2 (0 (3.15b) 

z 12 (i,r),(=0) = (1-1?) Z t (0 + t) Z 2 (0 (3.15c) 

The surface described by the above equations is known as a lofted surface between curves 1 
and 2 and is shown in figure 3-10(a). 

In a similar manner, we linearly interpolate between curves 3 and 4 to form a lofted 
surface described by 


(3.16a) 
(3.16b) 
(3.16c) 

This surface is shown in figure 3-10(b). 


*J < («.1.<=0) = (1-0 x 3 ( v ) + ( x 4 (0 
y M (0?,<= 0) = ( 1-0 y 3 (0 + ( Y 4 ( n ) 

*«({,>?, C=o) = (i-O z 3 W + < z 4 W 
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If we add these two surfaces, we obtain a surface which does not pass through the four 
corners determined by the intersections of the four twisted curves. In order to correct for this, 
we define a third surface by 


* 22 M«>1.<=°) = (l-O(l-i)) *,«=0) + { (1-!)) X 1 « = l) 

+ (1-0 n X 3 ((=0) + i i Jf,«=l) (317a) 

= (l— e)(i-v) y,u=o) + f (i-!») v,«=i) 

+ (1-0 I) Yg((= o) + ( n y a («=i) (3.17b) 

*1«<K.1,C= 0) = (1-0(1-!?) Zl((= 0) + f (1-!)) Z t ((= 1) 

+ (1-0 1? z 2 ((= 0) + 4 !) Zj(«=l) (3.17c) 

This surface interpolates the four corners of the region and is simply the plane section shown 
in figure 3-10(c). 

By adding equations (3.15) and (3.16) and subtracting equation (3.17) from the result, 
we obtain the following expressions which describe the surface bounded by curves 1 through 4 
(fig. 3-10(d)): 


*«.!>. c=0) = (1-!)) x,(0 + !? X 2 (t) + (1-0 * 2 (l) + i X 4 W 

- [(l-0(l->?) *,« = <>) + « (1-!)) ^l(f = l) 

+ (1-0 !) * 2 ({=0) + ( !) X 2 (« = l)] (3.18a) 

i«,1,C=0) = (1-1) Yj ( 0 + i) Y 2 (0 + (1-0 Y s (q) + { Y 4 ( V ) 

- [(i-o(i-i) y 2 (e= o) + « (i-!)) y 2 (f=i) 

+ (i-o i v 2 «=o) + « !) y 2 «=i)i (3.i8b) 

<?,!), <=o) = ( 1 - 1 ) z 2 (0 + 1 z 2 (0 + (1-0 z 2 (i) + ( ZM 

- [(l-O(i-i) Zj((= 0 ) + * ( 1 - 1 ) z,((= 1 ) 

+ (1-|) !) Z 2 «=0) + e 1 ^ 2 (« = l)l (3.18c) 
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Approximating a surface, by linear interpolation between curves that bound it, is 
known as transfinite bilinear interpolation. This method was named by Gordon and Hall (ref. 
34), but was first demonstrated by Coons (ref. 33). Thus, this method is also referred to as 
linear Coon’s interpolation. 

For some surfaces, transfinite bilinear interpolation does not produce a satisfactory 
approximation. This usually occurs when the curves that bound the surface possess extreme 
curvatures. In such cases, higher degree blending functions (e.g., cubic Hermite blending 
functions) can be employed and this is the subject of the next section. Also, if additional 
information about contour lines on the interior of the surface is available, one can break the 
surface up into simpler subsurfaces and approximate each subsurface separately. 


3.3.2 Three-Dimensional Bidirectional Hermite Interpolation. — As noted at the 
conclusion of the previous section, transfinite bilinear interpolation sometimes does not produce 
satisfactory surface approximations. Also, if a boundary surface is broken up into two or more 
subsurfaces with each subsurface generated by transfinite bilinear interpolation, then that 
boundary surface will have discontinuous first-order derivatives at all interfaces where different 
subsurfaces connect. Here, a new technique for generating surfaces, referred to as 3-D 
bidirectional Hermite interpolation, has been developed which does not have the 
aforementioned shortcomings of the transfinite bilinear interpolation. 

To illustrate the 3-D bidirectional Hermite interpolation, again consider the boundary 
surface bounded by four twisted curves shown in figure 3-8. Information about the four 
twisted curves may be given in analytical form, or they may be put into analytical form from 
coordinates of the nodal points located on the curves by using either cubic spline or tension 
spline interpolation described in Section 3.2. In either case, we end up with parametric 
equations describing each of the curves in the following form: 
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*,•=*<(«,•) Y i= Y i( s i) Zi=ZiM 


(3.19) 


where X '., Y t , and describe curve t and * = 1, 2, 3, 4. The parameter, s { , represents the 
approximate arc length along curve i as discussed in Section 3.3.1. 

Recall that for the Two-, Four-, and Six-Boundary Methods, boundary surfaces in the 
spatial domain are mapped onto coordinate planes in the transformed domain. Here, we map 
the boundary surface shown in figure 3-8 onto the coordinate plane shown in figure 3-9 located 
at £ = 0 in the transformed domain. Accordingly, we relate the s i parameters to £ and tj by 
equation (3.14). 

By using equation (3.14) and the Four-Boundary Method with Hermite interpolants 
described in Section 2.2, we obtain 

x(£,0, c=0) = X? 12 (t,r),(=0) + Az(£,ij,C=0) (3.20a) 

ytf,0,C=0) = i/' I2 ({,0>C=0) + Ay((,vX= 0 ) (3.20b) 

<f,i7,C=0) = ** I2 ((,nX=0) + Az«,q,C=0) (3.20c) 


where 


^(^,<=0) - x 2 (0 hfa) + X 2 (t) h 2 to) 

, <9x(£,7?=0,C=0) , , 5x(^,r;=l,<=0) , 

+ ^ h 3 (n) + ^ h 4 W 

;,C=0) = Ytf) hfa) + Y 2 (() h 2 (r)) 

. <9y(£,*?=0,C=0) L dytf, iy=l,C=0) , , * 

+ ^ ^ (T?) + ^ ^ 


(3.20d) 


(3.20e) 


* 12 {trt,C=o) = Z M) h M) + z 2 (0 h 2 ( v ) 

, dztf, 77=0,C=0) t dz(t,r,= 1,<=0) 

+ mi h3{r)) + 55 4{v) 


(3.20f) 
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Ax«,T),C=0) = Kf=0,>?,C=0) - x’j S ((=0,r/X=0)] h s (() 
+ [x(l = l,>),C=0) - i'j 2 (C=1,i;,C= 0)J h e (() 

. ,dx((=0,ri,(-0) 5x' J2 (l=0,v,C=0), , , c , 

+ l 51 M 1 

, ,<>*({ =m,c=o) a*' J2 (f=i,i?,c=o) 1 . 

+ 1 51 51 1 

A»(1,ij,C= 0) = (y(£=0,i),C=0) - y' I2 (?=0,y,C=0)J Ml) 

+ [ S /(€=l,r ? ,C=0) - s/' I2 (l=l,>?,C=0)] MO 

, ,5</(l=0,i?,C=0) 9y' 12 ((=0,riX=0). . . 

+ 1 51 5? 1 Mf) 

r 5y(l = l,rj,<=0) Sy' J2 (^ = l,i),C=0). . . . 

+ [ 1 MO 

A4t,n,(=0) = Kl=o, •»,<=(>) - *' j2 (€=0,ij,<=0)1 MO 

+ [«= M,<=0) - * 12 ((=l,t)X=0)] MO 

. ,d 2 ((= 0 ^x= 0 ) a* J2 (l=o, ?,(=(>), , 

+ 1 51 51 1 

i (52(1= l^C—O) 52' 22 (1=1,»7,1=0). . / t -. 

+ [ 5| 5^ ] MO 


ar' ; 2 (i=o,y,c=o) , e*(i= 0 ,i ? = 0 ,<= 0 ) , > 5 <i=o,y =i,c=Q) 

Sc = "ill) 51 + MW a/: 


51 


51 


, , ,a 2 x(l=0.')=0.C=0) , C 5 2 x(l = 0,v=l,<=0) 
+ MO aiay + MO a!5rj 


ax' 22 (1 — 1)^,1 — 0) A 5 x(l — 1 ,t; — 0,1 — 0) ^ ^ 5x(l Q) 

g^ MO 5^ + MO ai 


a 2 x(l=i.’?=0,<=0) 5 2 x(i=i ! 5=u=0) 

+ ^3\V) didr) a<*a„ 


didr) 


(3.20g) 


(3.20h) 


(3.20i) 


(3.20j) 


(3.20k) 
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dy'i2(£=°>»?>C= 0 ) 

di 


dy(t=0,T/=0,C=0) + h ^ <9y(^=0,r/=l,C=0) 


h M ^ 

*£ 

d^drj 


dt 


, A ,. '3*^=0, q=0,C=0) , h , > a p y(£=0,^=l,C=0) 

"b ' " 4 V u £»«?■£„ 


d^di) 


(3.201) 


^2/ , 2^(^ = 1 > r 7»C— °) i. /_x dy(£=l,»7=0,<=0) . r v 0y(£ = l,»7 = l,C=O) 

07 = M?) 37 + h 2iV) 37 


dt 


dt 


dt 


* ** r«~isr« * m .> 

»*'ji(f=o,i),<=o) _ , , a4«=o, 9 =o,c=o) , ,. , , d 2 (('=o,,=i,<=o) 

37 - 37 + n 2\V) 37 


^ - -IV'// ££ 1 ^ 

0*^=0, u = 0,C=0) . ,* ^^=0,^=1, C=0) 

+ h *w> a?^ + h4{v) Wv 


(3.20m) 


(3.20n) 


a* 12 (t= M,c=o) , , , a^=i,i7=o,c=o) , t &<e=M=u=o) 

^ = M*/) 37 + h 2\n) 37 


di 


dd 


9 2 ztf = l,r/=0,<=0) - . a 2 *({ = l ) ij=l,C=0) 

+ A *w 5 ^ + h “ M sm 


(3.20o) 


In the above equations, h 2 , h 2 , h 3 , and h 4 are given by equation (2.11) and h 5 , h ? , and 
are given by equation (2.24). The method for evaluating the partial derivative terms 
appearing on the right hand sides of equation (3.20) is described below. 

The first-order partial derivative terms in equation (3.20) determine the shape of the 
surface shown in figure 3-8 bounded by the four twisted curves given by equations (3.14) and 
(3.19). These derivative terms also determine whether grid lines on that surface will intersect 
the four twisted curves orthogonally or not. Since the only information we have about the 
surface shown in figure 3-8 is the four twisted curves, these curves are used to derive 
expressions for the first-order derivative terms in equation (3.20). 


We shall illustrate how first-order derivative terms in equation (3.20) are calculated by 
deriving expressions for and These three 

derivative terms are evaluated along curve 1 in which £ varies between 0 and 1, t)— 0, and C=0 
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(figs. 3-8 and 3-9). The method for deriving first-order derivative terms along the other three 
twisted curves (i.e., curves 2, 3, and 4 in fig. 3-8) will be similar to the procedure described 
below. 

Expressions for the first-order partial derivative terms along curve 1 are derived in five 
steps and they are 

1. Determine a vector which is perpendicular to the plane formed by the vectors 
tangent to curves 1 and 3 at one end of curve 1 (figs. 3-8 and 3-11). 

2. Determine a vector which is perpendicular to the plane formed by the vectors 
tangent to curves 1 and 4 at the other end of curve 1 (figs. 3-8 and 3-11). 

3. Linearly interpolate between the two vectors determined in Steps 1 and 2 

and denote this vector function as N^. 

4. Determine a vector function that is parallel to the cross product of (Step 

3) and a vector function tangent to curve 1. The vector function thus determined is 

assumed to be tangent to the surface that we wish to approximate along curve 1. 

5. Determine first-order partial derivative terms along curve 1 by forcing the 

vector function formed by the first-order derivatives to be parallel to the vector 
function determined in Step 4. 

The details of these five steps are described below. 

Step 1 - Determine Normal Vector at One End of Curve 1. The vectors tangent to 
curves 1 and 3 at £=0 and 77=0 are given by 

(3.21a) 
(3.21b) 

where I, J, and K are unit vectors pointing in the x-, t/-, and ^-directions, respectively. 


dX 2 ((=0) T , dY j(£=0) 1 , dZ^ = 0) 

dl I+ M J+ di 

dX 3 (r)=0) T , dY 3 ( V =0) T , dZ 3 (r)—0) 

~ foj 1 + dr) J + dr) 
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The vector perpendicular to the plane formed by the above two vectors is given by 


(13 X T„31 


' dY,((= 0) 

dz 3 ( n = 0 ) 

di 

dr) 

f dx,(i= 0 ) 

dZ 3 (r)—0) 

{ 

dri 

r dx t (i= 0 ) 

dY 3 ( n = 0) 

dt 

dr) 


az 1 («=o) 

dY 3 (i= 0) 

di 

dr) 

dZ 1 (i= 0) 

dX 3 ( V =0 ) 

di 

di] 

dY t (i= 0) 

o' 

II 

us 

>< 

dt 

dr] 


(3.22) 


Step 2 - Determine Normal Vector at Other End of Curve 1. The vectors tangent to 
curves 1 and 4 at £=1 and 77=0 are given by 


W^l) *Yi«=l) j , dZ,(t= 1) 

T ^4 “ — di 1 + dt. J + dt 

dx 4 (77=0) , dv 4 (77=0) T , az 4 (77=0) v 

T ^i = — ^ — 1 + d^q J + at K 


(3.23a) 

(3.23b) 


and the vector perpendicular to the plane formed by the above two vectors is given by 


X T r)41 


+ 


{ 

{ 

{ 


dY t (i= 1) 

dz 4 ( n = 0) 

dZ x (i = 1 ) 

Q? 

■erv 

II 

O 

di 

dr] 

di 

dt) j 

9X t (i= 1 

) dz 4 ( n = 0) 


9^4 (>7=0) 

di 

d<l 


dn 

dX 4 (i= 1 

) 9Y 4 (n=0) 

ay 2 u=:i ) 

dx 4 ((=0) 

dt 

dr] 


drj 


} 

} 


J 

K 


(3.24) 


Step 3 - Linearly Interpolate between Two Normal Vectors. In this step, we linearly 
interpolate between the two normal vectors obtained in Steps 1 and 2 to produce the following 
vector function: 
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N fi = (1 - 0 N 13 + Z N J4 


(3.25) 


where N i5 and N 14 are given by equations (3.22) and (3.24). 

Step 4 - Determine a Vector Function Tangent to the Surface. Since there are an 
infinite number of surfaces that can pass through any given curve, a strategy must be 
developed to determine which surface is to approximate the surface shown in figure 3-8. Here, 
the surface selected is the one that will pass through curve 1 as well as curves 3 and 4. Thus, 
the vector function tangent to the surface at curve 1 can be approximated by 

E„j = T fl X N {; (3.26) 


where N fi is given by equation (3.25) and T fi is a vector function tangent to curve 1 given by 




_ ggi(g r , dY M i x d JM 


dZ 


dz 


J + 


dZ 


K 


(3.27) 


Step 5 - Determine Expressions for First-Order Derivatives. Now that we know along 
which surface partial derivatives with respect to 7] at curve 1 can be made, we can derive 


expressions for the first-order partial derivative terms 

dz(Z,V= 0,c=0) 


dx(£, T?=0,C=0) dy(f,T/=0,C=0) 


, and 


dr] 


drj ’ dr] 

that appear in equation (3.20). Since we desire grid lines on the surface to be 


perpendicular to the boundary curve (curve 1 in this case), the vector function tangent to the 
^-direction, T ql , must be parallel to the vector function, E^, derived in Step 4. This can be 
expressed mathematically as 


T „ X E vl = 0 


(3.28) 


where 




dx(Z,v= 0,C=0) J dy(Z,ri= 0,C=0) j + dz(Z, t]=0X=0) r 


dr] 


dr] 


drj 


(3.29) 
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The desired expressions for the first-order partial derivative terms are obtained from 
equation (3.28) in a manner very similar to the derivation of equation (2.12). Again, “K 
factors” which appear in equation (2.12) can also be used here to ensure that grid lines do not 
overlap each other and to create a smoother surface. 


3.3.3 Parametric Bihyperbolic Spline Interpolation. — To illustrate how parametric 
bihyperbolic spline interpolation approximates surfaces, consider the organized set of IM x JM 
nodal points shown in figure 3-6. These points lie on surface H and are given on a rectangular 
xy-grid. We will call this collection of nodal points N and represent each nodel point within N 

by 


N(ij) ~ {Xy , Vy, 2y} , t = ^ ..., IM, j = ^ ..., JM 


(3.30) 


We wish to obtain a parametric description of surface H by interpolating between the 
nodal points of N. We will accomplish this in three steps. First, we will calculate approximate 
arc lengths along lines of constant i and j. Next, we will determine values of the derivatives at 
the nodal points which lie along the edges of surface H. Finally, we will use bihyperbolic spline 
interpolation to interpolate between the nodal points using our approximate arc lengths as the 
independent variables. This approach is similar to that suggested by Smith (ref. 36). 

We begin the first step by denoting approximate arc length along a line of constant i as 
s and approximate arc length along a line of constant j as t. We calculate s and t for each 
nodal point using the following expressions: 


s u — t n — 0 


(3.31) 


s i-i j + N 

1 (*«• 

t ij-i + \ 



~ x i-i j) 2 + (Vij-Vi-l j) 2 + i z ij~ z i-lj) 2 


x i j-l) “h (^tj j- j) “b ( z tj Z i j- 1 ) 


(3.32a) 

(3.32b) 
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where 


i = 2, 3, IM, j = 2, 3, JM 

Here, and t { - are the s and t values at nodal point N(i,j). 

Having accomplished this, the next step is to calculate values for the derivatives at the 
nodal points which lie along the edges of surface H. This information is required by the 
bihyperbolic spline interpolation method. The derivatives that are required are as follows: 

|f, If, and |f at i = 1, IM, j = 1, 2, JM (3.33a) 

If, If, and |f at N(iJ), i= 1,2,..., IM, j = 1, JM (3.33b) 

A & and A at jv(l ’ l) ’ w) ’ 

and N(IM,JM) (3.33c) 

We recommend using second-order-accurate finite- difference formulas to calculate these values 
from the coordinates of the nodal point coordinates. A comprehensive list of finite-difference 
approximations for derivatives appears in reference 60. 

Once we obtain values for the derivatives in equation (3.33), we are now ready to use 
bihyperbolic spline interpolation to form an approximation to surface H. Here, we follow 
closely Spath’s description of bihyperbolic spline interpolation found in reference 50. In 
bihyperbolic spline interpolation, we must determine the coefficients a { - kl of the expression 

u(s,t) = X/ S a ijkl ^k( S ^ a i^ S ij' S i+l j) j) (3.34) 

k=l 1=1 

where u represents either x , y, or z in the Cartesian system. The function u(s,t) in equation 
(3.34) must take on the specified values of x , y, or z (depending on which one it represents) at 
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the appropriate nodal point locations while maintaining continuity 

up to the second-order 

derivative everywhere. The functions <j> k and <j>i are given by 


<*,•» s ij > s i+i j) = s - s ij 

(3.35a) 

a i> s ij-> s i+i j) ~ s i+l j~ s * 

(3.356) 

s ij' s i+l j ) = V’( 5— s ij’ a i' S ij » S t l+l j ) 

(3.35c) 

a »» S ij' S i+1 j ) = S i+1 j~ S ' a i ’ S ij' S i+1 j ) 

(3.35d) 

where 


A s..stnft(«.i) — s stnh(ajAs--) 
S ’ *«’ ®*'+J ^ _ smA(a t As 0 ) - 

(3.35e) 

As^- = s i+i j~ s ij 

(3.35f) 


Here, a t - and /?■ are tension parameters in the s and t “directions,” respectively. They are 
constrained by 0 < < oo and 0 < /?■ < oo with tension increasing as they approach 

infinity. Each grid cell (fig. 3-6) can have different oc i and /? • values depending on the amount 
of tension desired. Here, we note that as a ■ and /?• approach zero, the bihyperbolic spline 
surface approaches a bicubic spline surface (ref. 50). 

The reader may recognize that equation (3.35e) contains expressions similar to those in 
equation (3.12) of Section 3.2.2. In fact, equation (3.12) can be cast in the form 


X i( s ) = Y A ik <t>k( s '<*i' s i' s i+i) 
k= 1 


(3.36) 


where <f> k is given by equation (3.35), and the A ik s are the coefficients which contain the 
unknown X"(s t ) values. Thus, the bihyperbolic spline described by equation (3.34) is merely 
an extension of the tension spline of Section 3.2.2. 


We begin the algorithm for obtaining the coefficients, a - ki , by letting u - = u(s-, t-), 
du • du • d 2u ij 

p (J = = -{ft' an< ^ r ij = q s " Qf R- ecaI1 that u is known at all nodal points, while p-, 

q- , and r- are only known at the nodal points which lie on the boundaries of H. Values of p 
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q and r- at nodal points other than those on the boundaries can be determined by solving 
the following 2(/M) + JM + 2 linear system of equations: 


Here, v- 


^i-1 j Pi-1 j + j V i-1 j + Pij + j 

l ij~ U i-l j 




= Ku K-o + 

for j = 1, 2, . . JM , i = 2, 3, . . IM—1 


c i j- 1 Qi j-1 ( C t j- 1 W i j-1 c ij W ij) Vij + C ij j+1 

“ H+l-% 


= <ii-i to* + + ^^str 


-» j-i 

for i = 1, 2, /M , j — 2, 3, JM— 1 


(3.37a) 


(3.37b) 


[ j r t'-J j 

+ (Vi 

■ u- + £>•• u 

3 i-l 3 *3 i 

= W.J 

i ( U .-I 3 

^ S t-2 j 

for j = 

1, JM , 

i — 2, 3, 

-1 r ij-l 

+ 

i “ij-I + C H *. 

= C ij- 

2 K>J 

+ 

^ l i 3-1 


v^t j+1 Pij 


for i — 1, 2, ..., /M , j = 2, 3, JM— 1 
j, w ijt b-, and c {j are given by 


tv- 


j+l 


6 ,, = 


a • 2 A $ - smh( a ■ A s^- ) 


U (%• 2 -l)(smh(a | As ij )-a | As. i ) 


C • • = 


(3 j 2 &tjj sin h((3 iAtjj) 


tJ (w^-XKsinhiPiAiJ-piAtJ 


(3.37c) 


JM— 1 


(3.37d) 

/• = 1, 2, . 

IM—l) 

(3.38a) 

f 0‘= 1.2,-. 

JM—1) 

(3.38b) 

(»’ = 1, 2, . 

IM—1) 

(3.38c) 

0 = 1, 2, . 


(3.38d) 
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The system of equations described by equations (3.37) and (3.38) is obtained by requiring that 
the function u and its partial derivatives p, q, and r be continuous at all of the interior nodal 
points in N. This system has a tridiagonal coefficient matrix that can be shown to be 
diagonally dominant (ref. 50). Thus, the system may be easily solved for the unknown values 
of p-, q-, and r- by using the Thomas algorithm. 

We now solve for the a^’s of equation (3.34) by first defining a 4x4 coefficient matrix 
within each grid cell as 


A; j — Uijki (k, l — 1, 2, 3, 4) 

This matrix can be found from the following matrix product expression (ref. 50): 
The variables C and K { - are matrix expressions which are given by the following: 


(3.39) 


(3.40a) 


C(g,h) = 


K ij = 


0 

0 

1/^ 

0 

1 /g 

0 

0 

0 

l 

-1 

1 

-h 

g( l-h) 

t- 1 

1 

5r 

to 

g( i-fc) 

1-h 2 

1 

h 

l 

1 

g(l-h) 

1 -h 2 

S' 

i 

ji-H 

i 

1 -/ 1 2 

u ij 

q ij 

U i 3+1 

q i j+1 

Prj 

r ij 

P ij+1 

r ij+l 

U i + 1 j 

q i+l 3 

U i+1 3+1 

q i+l j+1 

P i+1 j 

r i+l 3 

P i+1 j+1 

r i+l j+1 


(3.40b) 


(3.40c) 
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Having found the coefficients we can now use equation (3.35) to calculate w(s,<) at any s 

and t. 

Returning to the problem of obtaining a description of surface i/, we perform three 
bihyperbolic interpolations (one for u=x{s^t), one for u—y(s,t), and one for u—z(s,t)) using the 
procedure above. In this way, we obtain the parametric representation of surface H that is 
necessary for the grid generation technique presented in Section 2.0. 
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4.0 GENERATION OF SINGLE AND COMPOSITE GRIDS 


In the previous two sections, the details of the algebraic grid generation methods used 
in GRID2D/3D are described. In this section, we demonstrate the usefulness of GRID2D/3D 
by using it to generate a number of single and composite grids within complex-shaped 2- and 
3-D spatial domains. Recall from Section 1.0, a single grid is a grid system based on one 
boundary-fitted coordinate system, and a composite grid is a grid system made up of two or 
more single grids patched together. 


4.1 Single Grids 


All of the details involved in using the Two-, Four-, and Six-Boundary Methods to 
generate single grids were presented in Sections 2.0 and 3.0. In fact, while illustrating these grid 
generation methods, several single grids were generated (figs. 2-3, 2-4, and 2-8). In this 
subsection, several additional single grids generated by GRID2D/3D are presented to illustrate 
the capabilities of GRID2D/3D. Also presented in this section is a brief discussion on how to 
smooth discontinuities in grid systems that arise from boundary discontinuities. 

Figure 4-1 shows a 2-D single grid generated for an irregular, 2-D, coastline topology 
using the Two-Boundary Method. A stretching function was used to cluster grid points near 
the coastline where complex flow conditions are expected. The coastline curve was formed by 
using parametric tension spline interpolation. 

Figure 4-2 shows a 2-D single grid around a sharp bend generated by using the Two- 
Boundary Method. This figure illustrates how boundary discontinuities can propagate into the 
interior of the grid. The slope discontinuity of the r) grid lines along i = lj can be eliminated 
by applying the following equations: 
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(4.1a) 


v “ r = yh + - 25 M + U - 2 + y”-ij) 


<r = <J + .25 (*f +1|i - 2 xlj + *SL lpi ) 


(4.1b) 


where i = I 2 — (m can be zero or some positive integer constant) and j = 

2,3,...,JL — 1. In order to smooth the grid shown in figure 4-2, equation (4.1) needs to be 
applied a number of times; that is, n = 0, 1, 2, 3, ... with n = 0 being the grid shown in figure 
4-2, n = 1 being the first correction, n = 2 being the second correction, and so on. By using 
equation (4.1) repeatedly, the grid shown in figure 4-2 was smoothed as shown in figure 4-3. 

Figure 4-4 shows a 2-D single grid around an airfoil. This grid was generated by using 
the Two- Boundary Method in which stretching functions were used to cluster grid points near 
the airfoil surface in anticipation of the complex boundary layer flow there. 

Figure 4-5 shows a 2-D single grid for an irregularly shaped, four-sided region. That 
figure clearly shows how orthogonality of the grid lines is enforced at each of the four 
boundaries of the spatial domain. This grid was generated by using the Four-Boundary 
Method and stretching functions were not used. The four boundary curves of the spatial 
domain shown in figure 4-5 were represented in parametric form by tension spline 
interpolation. 

A 3-D single grid between the twisted blades of a radial turbine is shown in figure 4-6. 
The Two- Boundary Method was used to generate this grid and a stretching function was used 
to cluster grid points near the blade surfaces. The blade surfaces were generated by using 
tension splines in conjunction with bilinearly blended transfinite interpolation. 

As an example of the efficiency of GRID2D/3D, the typical real time taken to generate 
a 21 x 21 x 21 3-D single grid is 45 seconds on an IBM AT compatible personal computer 
running at 12 MHz. 
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4.2 Composite Grids 


In Section 1.1, it was noted that GRID2D/3D can be used to generate composite grids 
as well as single grids. Also, it was noted that composite grids which can be generated by 
GRID2D/3D include those that are completely discontinuous, partially discontinuous, and 
partially continuous (fig. 1-2). For completely or partially discontinuous composite grids such 
as those shown in figures. l-2(a) and l-2(b), each single grid within the composite grid can be 
different from any other in structure and in the number of grid points. Thus, each single grid 
within such composite grids can be generated independently of the other single grids, and each 
single grid can be generated in the manner described in Sections 2.0 and 4.1. How the 
different single grids of such composite grids are patched together once they are generated is 
described in Sections 1.1.1 and 1.1.2. Since discontinuous composite grids can be generated 
rather easily and patching is trivial, no further discussions concerning them will be given. 

For partially continuous composite (PCC) grids, each single grid within it must satisfy 
a number of compatibility conditions so that when the different single grids are patched 
together, the resultant composite grid will have some degree of continuity (Section 1.1.3). 
PCC grids such as the one shown in figure l-2(d) are often difficult to generate. However, 
when it is possible to generate them, they are the easiest to use with FD and FV methods to 
obtain solutions to partial differential equations. Below, we discuss how PCC grids can be 
generated. 

4.3 Partially Continuous Composite Grids 

PCC grids are generated by GRID2D/3D in two major steps. The first step involves 
partitioning the spatial domain into zones. The second step involves constructing grid systems 
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that will be continuous from one zone to another. These two steps are described in detail 
below. 


4.3.1 Partitioning. — The first step in constructing a PCC grid is to partition the 
spatial domain of interest into a finite number of contiguous zones. The partitioning process 
involves answering three interrelated questions: 

1. How should the spatial domain be partitioned into zones? 

2. Based on that partition, what grid structure is to be used within each zone? 

3. Based on that partition and grid structure selections, how should the 
different zones be mapped to the transformed domain so that the resultant 
composite grid will have some degree of continuity? 

The answers to these questions depend on the geometry of the spatial domain and the 
physics of the problem for which the grid is being generated. For any given problem, several 
different choices are usually possible. However, for simplicity, the choice containing the least 
number of zones is often the most desirable. 

An example illustrating one strategy at partitioning is shown in figure l-3(a). The 
boundaries of that spatial domain contain a backward facing step, a flat wall, and a wedge- 
shaped obstacle. The fluid is flowing from left to right. In figure 1-3, the spatial domain is 
partitioned by using free-streamline theory. More specifically, a new zone is set up where 
separation is expected (e.g., at the backward step and at the rear of the wedge-shaped 
obstacle) or where a flow will separate into two streams (e.g., at the nose of the wedge-shaped 
obstacle). The structure of the grid within each zone in figure l-3(a) was chosen to be H-type. 
This structure aligns the main flow direction with the grid lines which is important because it 
reduces dissipation error and permits the use of the thin-layer Navier-Stokes equations. Also, 
this structure allows grid lines to be clustered near solid wall boundaries readily as well as 
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allows continuity of grid lines and their derivatives from one zone to another. Figure l-3(b) 
shows how the different zones in figure l-3(a) are mapped to the transformed domain. 

As a second example, consider the 2-D inlet of an aircraft engine shown in figure 4-7. 
This geometry would be awkward to handle with a single grid due to the separation of the 
internal flow from the external flow by the cowl. Figure 4-8 shows a partitioning of this 
spatial domain into two zones and the transformed domain to which it is mapped. Other 
configurations involving more than two zones are certainly possible, but more complicated 
configurations are unnecessary unless the physics of the flow dictates their use. If we knew 
something about the flowfield within this spatial domain, we might wish to introduce new 
zones in order to better resolve the physical aspects of the flow (e.g., to make capturing of 
shock waves more convenient). In some cases, we might even want to solve different 
governing equations in different zones in order to obtain solutions more rapidly. 

As a third example, consider the spatial domain shown in figure 4-9. This is a 2-D 
region near a wind tunnel wall with slots cut into it. Two possible partitions for this spatial 
domain are shown in figure 4-10. The partition shown in figure 4-10(a) contains six zones. A 
better partition with only three zones is shown in figure 4- 10(b). 


4.3.2 Grid Generation with Patching. — Once we have partitioned the spatial domain 
of interest into zones, we are ready to generate a grid for each zone. For partially continuous 
composite (PCC) grids, the grid generated for each zone must be such that when they are 
patched together, the resultant composite grid has some degree of continuity. It turns out 
that this requirement necessitates some minor modifications to the Two-, Four-, and Six- 
Boundary Methods presented in Section 2.0. In this section, these modifications are described 
in the framework of the Two-Boundary Method by applying it to generate a PCC grid in the 
2-D inlet of the aircraft engine shown in figure 4-7. 


66 



Step 1 - Define the Coordinate Transformation. We seek a coordinate 

transformation of the form given by equation (2.19). 

Step 2 - Select a Time-Stretching Function. We set r equal to t according to 
equation (2.2). 

Step 3 - Select Two Boundaries of Each Zone. We choose curves 1 and 2 for 
zone 1 and curves 5 and 6 for zone 2 (fig. 4-8). For each zone, we map these two curves to 
coordinate lines V = V low and V = n high , respectively. Here, rj low , r) high , and £ high are 

defined as the coordinate lines in the transformed domain which correspond to the boundaries 
of a zone. Thus, for zone 1 we have 


x, = *«,<>=>), 0 j = *j«) 

y, = *(,*> =v h j = yj«) 
x 2 = <t*=ri hiah ) = -Vi(f) 

y 2 = viZ,r,= nhigh ) = Y 2 (i) 


(4.2a) 

(4.2b) 

(4.2c) 

(4.2d) 


where rj low = 0 and r) higfl = 0.5. and Y 1 are the x- and y-coordinates of curve 1, and X 2 and 
Y 2 are the x- and y-coordinates of curve 2. 

For zone 2, we have 


*5 = = X S (Z) (4.3a) 

Ys = ! = = Y s (() (4.3b) 

x e = x ti>'l = ’> Mgll ) = x eti) (4.3c) 

Ye = vti,<l=n high ) = Yeti) (4.3d) 


where r) [ow = 0.5 and r) hih =- 1. X 5 and Y 5 are the x- and y- coordinates of curve 5, and X 6 
and Y 6 are the x- and y-coordinates of curve 6. 
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The other two curves in each zone (curves 3 and 4 in zone 1, and curves 7 and 8 in 
zone 2) are mapped to coordinate lines £ = t low and t = thigh in the transformed domain. For 
our example, ti ow = 0 and thigh.— ^ ^ or both zones 1 and 2. 

If we were using the Four- or the Six- Boundary Method, then we would have selected 
four or six boundaries in each zone as described in Sections 2.2 and 2.3. 

Step 4 - Describe the Two Boundaries of Each Zone in Parametric Form. We now 
represent the four curves selected in Step 3 (two for each zone) in parametric form. Here, 
information about the four curves is given in the form of four sets of discrete points which lie 
along the curves. Thus, we use the 2-D parametric tension spline interpolation method 
described in Section 3.2.2 to generate the parametric equations in the form dictated by 
equations (4.2) and (4.3). By using this method, we generate the 2-D curve approximations 
shown in figure 4-11. 

Step 5 - Define Curves That Connect the Two Boundaries in Each Zone. As before, we 
use transfinite interpolation to derive curves which connect the two boundaries of each zone 
described in Step 3. For Lagrange interpolation, the equations for the connecting curves 
between two boundaries in any one zone are 


4t>v) = x (trt = Vi ow ) h^) + x it,v=Vhi 9 h) hi 1 )) 
y(t,T}) = y(t,v=Vi 0W ) hiv) + y(t^=Vhigh) hi 1 !) 


where the blending functions, and l 2 , are given by 


l liv) 


V-Vhigh 

*7 low *7 high 


hi 1 !) = 


V-Vlow 

*7 high *! low 


(4.4a) 

(4.4b) 


(4.4c) 


(4.4d) 
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As stated in Section 2.1, the curves described by equation (4.4) are straight lines which 
do not, in general, intersect the boundary curves perpendicularly. This could lead to slope 
discontinuities in grid lines which cross zonal interfaces (fig. 4-12). For this reason, we 
recommend the use of transfinite interpolation based on Hermite interpolation when patching. 
Hermite interpolation allows the connecting curves to intersect the boundary curves 
perpendicularly. Thus, grid lines in adjacent zones which meet at the zonal interface will 
possess identical slopes at the interface. This eliminates slope discontinuities of the grid lines 
at zonal interfaces. 


For Hermite interpolation, the functional form for the connecting curves between two 
selected boundaries of a zone is as follows: 


4Z,V) = <t,V = how) M 7 ?) + x(t,V -Vhigh) h 2(v) 


+ g f tt ! *p ML ) hM + 


dxj^rj = rj low ) 
dr) 


h 4 (v) 


y(Z,v) = y(^v=viou,) h i(v) + y(^v=Vhi 9 h) M 7 ?) 


+ 


hm 


In the above equations, h 2 , h 3 , and h 4 are constrained by 


(4.5a) 


(4.5b) 


h l(v = View) = 1 

h l(V = Vhigh) = 

dh l(V = 1 UoJ_ n 
drj u 

dh l(V = Vhigh) _ 
dr) 

M 7 ? = v lo J = 0 

M 7 ? = Vhigh) = 

dh 2 (V = Vlow) _ n 
dr) 

dh 2 (r) = r) high ) _ 
dr) 

h 3 (V = Vlow) = 0 

h 3 (V = Vhigh ) = 

dh 3 (V = V io J _ 
dr) 

dh s(v = Vhigh) . 
dr) 
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h 4 ( r l = *Uow) = °» 


h 4 (*? = = 0 


= ’Jw) _ „ = ^-gjj _ , 

55 “ a? 

With these constraints, /i h 2 , and h 4 become 

(4.6a) 
(4.6b) 
(4.6c) 
(4.6d) 

where the following expressions hold: 


h 1 = atf 3 + btf 2 + ctf + dj 

h 2 = + b 2 rj 2 + c 2 t? + d 2 

h 3 = a 3 r) 3 + b 3 r) 2 + c 3 r) + d 3 

h 4 = + 6 4 t/ 2 + c 4 r) + d 4 


2 

Si*] low "Vhigh)^ low high) + low*] high. — *^ 1 ! low ~ r l high) 


(4.6e) 


l(V low *lhigh) 
2 (*1 low *1 high) 


(4.6f) 


°i = ~3*1 iow a i - 2r how b i 

d l — — *lhigh a l — *lhigh b l ~ *lhigh°l 

a 2 ~ ~ a l 



d 2 ~ 1 — d-L 


(4.6g) 

(4.6h) 

(4.61) 
(4-6j) 
(4.6k) 

(4.61) 


(*1 low *1 high) 


^(vfow — ^high^low *1 high) + 2(3rj low*! high 2 *llow *lhigh) 


(4.6m) 


b _ 1 ~ 3a 3 (*1 low -*l high) 

2 (*llow *1 high) 

c 3 = 1 — 3r) 2 QW a 3 — 2 Vi ow b 3 

d 3 ~ ~ *lhigh a 3 ~ *lhigh b 3 ~ *1 high 0 3 


(4.6n) 


(4.6o) 

(4.6p) 


a 4 — a 3 


(4.6q) 
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(4.6r) 


b _ I " 3 ^(vfow-Vhigh) 

low high) 

c 4 = 1 - 3 rjf ow a 4 - 2r} hw b 4 (4.6s) 

d 4 = ~Vhigh a 4 - r lhigh h 4 ~ 'l high 0 4 ( 4 ‘ 6t ) 

Similar to the examples in Section 2.0, the partial derivative terms in equation (4.5) are 
chosen so that the connecting curves will intersect the two boundary curves perpendicularly. 
Thus, for zone 1 we have 


dx{Z,T1 = T1 low ) _ 
dy 


(4.7a) 

dy(^V=Vi 0 J _ 

dr] 


(4.7b) 

9z(t,V = V h igh) _ 
dy 


(4.7c) 

dy(Z,n=ii l igh) _ 

dr) 


(4.7d) 


where r) low = 0 and rj higfl — 0.5. Likewise, for zone 2 we have 


dx(t,r) = r) low ) 
dr) 


(4.8a) 

dy(Z>v=Vi ow ) _ 
dr) 

" a f 

(4.8b) 

dx(Z,t)z=r) high ) _ 
dr) 


(4.8c) 

dy(Z,v=y high ) _ 

dr) 


(4.8d) 


where r) low = 0.5 and t) high = 1.0. 

An alternative procedure to the one just described for calculating and y(f,rj) is 

to do two mappings for each zone. The first mapping is from (x,y) to (^ , ,^/ , ) with boundaries 
of both C and y' between 0 and 1. The second mapping is from (£\r/’) to (£,t?) with the 
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boundaries of £ between £—Z low and £=f hi h and the boundaries of r) between V~Vi ow anc ^ 
Tj=Tjf l - gh . The first mapping can be accomplished by using the equations in Section 2.0 with £ 
and t) replaced by C and rj\ The second mapping is given by £ = £ low -f £’ H^igh ~ £ low ) 
an d 1 = V, ow + n’ (.V h , gh ~ r, low ). 

Step 6 - Discretize the Domain. We discretize the domain in the £- tj-t coordinate 
system by replacing the temporal domain with equally incremented time levels and by 
replacing each zone of the spatial domain with equally spaced grid points (fig. 4-13(b)). The 
time levels are described by equation (2.14). For our two-zone example, the grid points in 
zone 1 are located at (£,-, 77 ^), where 


£ t - = («- l)Af, t=l, 2, .... IL, (4.9a) 

r)j = 0—1) A i? , j=l, 2, JLl, (4.9b) 

and 

A ^ = 7lb 

The grid points in zone 2 are located at (£,-,7?j), where 

(i = (i-l) A? , «= 1, 2, IL (4.9d) 

r)j = 0—1) A , j= JLl, JL1 + 1,..., JL (4.9e) 

By substituting equation (4.9) into equation (4.5), we obtain the locations of the grid 
points in the x-y-t coordinate system (fig. 4-13(a)). 

We note that grid points IL1 through IL along grid line JLl in the transformed 

domain represent two grid points in the spatial domain. When using this grid to obtain 

numerical solutions, special care must be exercised to ensure that the correct grid point is used 
at the appropriate situation. 


A ij = 


JL - 1 


(4.9c) 
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Step 7 - Control the Distribution of Grid Points in Each Zone. For high-speed flows 
through the inlet, we expect large velocity gradients next to all solid surfaces. Thus, grid 
points need to be clustered near solid surfaces. With this in mind, we use stretching functions 
to cluster grid points near curves 1 and 2 for zone 1 and near curve 5 for zone 2. The new 
distribution of grid points is shown in figure 4-14. 

Here, we make a few comments about the continuity of grid lines and grid spacings 
across zonal interfaces. In the above example, grid lines remained continuous along the 
portion of the zonal interface which does not represent a “physical” boundary. This is a 
desirable situation since it simplifies the numerical algorithm which will be used to investigate 
the flow in the spatial domain. Proper allignment of grid lines at the zonal interface in this 
case is made possible by using an appropriate stretching function in each zone. Thus, one 
should exercise care when constructing a grid to ensure that grid lines remain continuous 
across the sections of the zonal interfaces which touch each other in the spatial domain. For 
our example, we also note that the grid spacing in the rj-direction must not change too 
abruptly across the zonal interface. Grid spacing in the ^-direction is defined as the distance 
between adjacent grid points along a grid line of constant As with grid alignment, proper 
grid spacing at the zonal interface depends upon using an appropriate stretching function in 
each zone. 

Step 8 - Calculate Metric Coefficients. For the spatial domain shown in figure 4-7, 
the five metric coefficients which need to be evaluated are r t , £ x , £y an d *?y These metric 
coefficients can be determined by using equation (2.28). We note that one-sided differencing 
should be used when calculating partial derivative terms at grid points which lie along the 
portions of the zonal interface which do not touch in the spatial domain. Central differencing 
can be used at the other grid point locations along the zonal interface. 
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5.0 SUMMARY 


A computer program — referred to as GRID2D/3D has been developed to generate 
single and composite grid systems in complex-shaped 2- and 3-D spatial domains. This 
technical memorandum describes the details of the algebraic grid generation methods used in 
GRID2D/3D, namely, the Two-, Four-, and Six-Boundary Methods. This technical 
memorandum also describes the methods used to derive parametric representations of curves 
and surfaces needed by the grid generation techniques. Several grid systems generated by 
GRID2D/3D for various 2- and 3-D spatial domains were presented to illustrate that 
GRID2D/3D is efficient and capable of handling complex geometries. 
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(A) 



<B) 



(C) 


(a) STRUCTURED. 

(b) UNSTRUCTURED. 

(c) MIXED (PARTLY STRUCTURED, PARTLY UNSTRUCTURED). 
FIGURE 1-1. - TYPES OF GRID SYSTEMS. 
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(D) 

(a) COMPLETELY DISCONTINUOUS. 

(b) PARTIALLY DISCONTINUOUS. 

(c) PARTIALLY CONTINUOUS. 

(d) PARTIALLY OR COMPLETELY CONTINUOUS. 
FIGURE 1-2. - TYPES OF COMPOSITE GRIDS. 
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(a) IN SPATIAL DOMAIN (x-y COORDINATE SYSTEM). 

(b) IN TRANSFORMED DOMAIN (£-f] COORDINATE SYSTEM). 

FIGURE 1-3. - CONTINUOUS COMPOSITE GRID SYSTEM MADE UP OF FOUR SINGLE GRIDS. 



FIGURE 1-4. - PARTIALLY DISCONTINUOUS COMPOSITE GRID THAT APPEARS LIKE A PARTIALLY OR 
COMPLETELY CONTINUOUS COMPOSITE GRID. 
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(a) 


(b) 


(a) IN SPATIAL DOMAIN, 
(b) IN TRANSFORMED DOMAIN. 

FIGURE 2-3. - GRID SYSTEMS. 



F IGURE ?-'l. - GRID SYSTEM IN THE SPATIAL DOMAIN AFTER STRETCHING. GRID IS COM 
PIJTER GENERATED. 
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(b) 


(a) IN x-y COORDINATE SYSTEM. 

(b) IN COORDINATE SYSTEM. 

FIGURE 2-5. - SPATIAL DOMAIN. 



FIGURE 2 6. - APPROXIMATIONS OE CURVES 1, 7. 3. AND <T OBTAINED BY USING 
PARAMETRIC TENSION SPLINE INTERPOLATIONS. THE CURVES ARE COMPUTER 
GENT RATED AND NODAL POINTS ARE DENOTED BY 0. 
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FIGURE 2-7. - MAPPING BE I WEEN THE SPATIAL AND TRANSFORMED DOMAINS OBTAINED BY USING EQUATION ( 2 . 21 ). NOTE THAT CURVES 3' AND 'J' IN THE SPATIAL 
DOMAIN ARE MAPPED TO COORDINATE l INES { - 0 AND $ = 1 IN THE TRANSFORMED DOMAIN, RESPECTIVELY. 





Spline Parameter, s 

FIGURE 3-1. - NOTATION USED IN DERIVING PARAMETRIC CUBIC SPLINES. NODAL POINTS ARE DE- 
NOTED BY o. 




(a) CYCLIC END CONDITIONS. 

<b) NATURAL END CONDITIONS. 

EIGURr 3?. - EXAMPLES OE CURVES GENERATED BY USING PARAMETRIC CUBIC SPLINE INTERP0LA1 ION. THE CURVES ARE COMPUIER GEN- 
FRAM D AND NODAL POINTS ARF DENOTED BY ►. 


OF 
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(a) CURVE BEING APPROXIMATED, 
(b) CUBIC SPLINE APPROXIMATION. 


FIGURE 3-3. - CASE WHERE A CUBIC SPLINE POSSESSES UNDESIRABLE WIGGLES. CURVE IN (b) IS COMPUTER GENERATED AND NODAL POINTS ARE 
DENOTED BY ► . 




C.O t , , | i ] . , . 1 

0*2 :.C 2.0 3.0 ft.O 5.0 

FIGURE 3-U. - A PARAMETRIC TENSION SPLINE WITH 0 = 10 FOR THE 
SAME SET OF NODAL POINTS AS THE CUBIC SPLINE OF FIGURE 3-3(b> 
NOTE THAT THE WIGGLES PRESENT IN THE CUBIC SPLINE SHOWN IN 
FIGURE 3-3(b) HAVE BEEN REMOVED. THE CURVE IS COMPUTER GEN- 
ERATED AND NODAL POINTS ARE DENOTED BY ►. 
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FIGURE 5 5. - EXAMPLES OF SURFACES WHOSE INTERIOR REGIONS ARE "SIMILAR" TO THEIR BOUND- 
ARY CURVES. BOUNDARY CURVES ARE SHOWN AS SOLID LINES, WHILE INTERIOR CURVES ARE SHOWN 
AS DASHED LINES. 



FIGURE 36. - POINTS ON SURFACE II GIVEN ON A RECTANGULAR X-y GRID. NODAL POINTS ARE 
DENOTED BY o. 





FIGURE 3-7. - EXAMPLES OF SURFACES WHOSE INTERIOR REGIONS DIFFER GREATLY FROM THEIR BOUNDARY CURVES, 


BOUNDARY CURVES ARE SHOWN AS SOLID LINES, WHILE INTERIOR CURVES ARE SHOWN AS DASHED LINES. 



FIGURE 38. - FOUR TWISTED CURVES WHICH INTERSECT TO FORM A CURVILINEAR "QUADRILATERAL." 
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FIGURE 4-2. GRID SYSTEM GENERATED ABOVE A SHARP BEND BY USING THE TWO-BOUNDARY METHOD. GRID IS COMPUTER GENERATED. 



FIGURE 4-3. - GRID SYSTEM OF FIGURE 4-2 AFTER SMOOTHING. GRID IS COMPUTER GENERATED. 
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FIGURE '15. - GRID GENERATED EOR AN IRREGULARLY SHAPED, 
GFNf RAIED. 



OUR-SIDED REGION BY USING THE FOUR-BOUNDARY METHOD. GRID IS COMPUTER 
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(a) IN x-y COORDINATE SYSTEM. 

(b) IN £-n COORDINATE SYSTEM. 

FIGURE '1-8. - PARTITIONING THE SPATIAL DOMAIN INTO /ONES. 






y 




(a) PARTITION WITH SIX ZONES. 

(b) PARTITION WITH THREE ZONES. 

FIGURE <4-10. - TWO POSSIBLE PARTITIONS FOR THE SPATIAL DOMAIN SHOWN IN FIGURE 4-9. 



FIGURE 4-11. - APPROXIMATIONS OF CURVES 1. 2, 5 AND 6 OBTAINED BY USING PARAMETRIC 
TENSION SPLINE INTERPOLATION. THE CURVES ARE COMPUTER GENERATED AND NODAL POINTS 
ARE DENOTED BY o. 
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boundary curve a 



rir.URF 'I 12. - AN FXAMPLF OF Sl.OfT D T SCON f I Nil ITIF.S PRESENT IN GRID LINES 
WHICH CROSS A 70NAI IN IF REACT . 



(a) IN THE SPATIAL DOMAIN, 
(b) IN THE TRANSFORMED DOMAIN. 


FIGURE R-13. - GRID SYSTEM. 






FIGURE '1-14. - GRID SYSTEM IN THE SPATIAL DOMAIN AFTER STRETCHING. GRID IS COMPUTER GENERATED. 
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